MVA Course
###install required packages###
#delete hash symbol on the following lines and run (only needs to be done first time, not every time you run the code afterwards!)
#install.packages("vegan")
#install.packages("dplyr")
#install.packages("corrplot")
#install.packages("cocorresp")
#install.packages("tibble")
#install.packages("labdsv")
#install.packages("knitr")
#install.packages("dendextend")
#install.packages("pls")
#install.packages("adespatial")Exercise 1. Introduction to R Notebooks and Markdown.
This exercise is only to get used to working with R Notebook and R Markdown. The data in this first are not multivariate and the statistical analysis is a simple linear regression. All following exercises will start with this kind of introduction.
For each exercise there will be a number of explicit questions written in green.
To pass the course, you should give your answers to these questions to one of the teachers during the computer labs.
Write your answers to the questions and any other comments here below.
Replace this text in for each question with your answers and comments…
A quick note on R Notebooks and R Markdown
These terms are often used interchangeably to mean a combination of R code, results, and formattable text, such as this one. They allow you to present your results, the code used to generate the results and your interpretation and comments all together in one document that you can easily update, and which is easy for a reader to understand (R Notebooks are a minor update of R Markdown files that allow for easier editing as you can run just part of your code and see the results in a preview. Both notebook and markdown are .Rmd files, and are written in exactly the same way, so don’t worry about the difference for this course!).
Example of a simple analysis
First we load or create the data we will be using. The box below is an example of a code chunk.
df <- read.table(header = TRUE, sep = ",",
text = "x,y
1,2
1.5,3.5
2,5
2.5,5
3,7")
# Anything you write in a code chunk after a hash mark becomes a comment like this
# and is ignored when running the code. If you are working in an R markdown notebook
# such as this however, you can also write your longer notes in the text sections
# outside the code chunks.So I could write my notes here instead, which is probably easier to read for anything more than a few words of explanation.
Now we can do a simple linear regression. The results are stored in the object “reg”.
reg <- lm(y~x, df) Plot the result.
plot(y~x, df)
lines(df$x, predict(reg), col = 'red') # Add a regression line to the plot
Get a summary of the results
summary(reg)##
## Call:
## lm(formula = y ~ x, data = df)
##
## Residuals:
## 1 2 3 4 5
## -0.20 0.15 0.50 -0.65 0.20
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1000 0.6819 -0.147 0.89271
## x 2.3000 0.3215 7.155 0.00562 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5083 on 3 degrees of freedom
## Multiple R-squared: 0.9446, Adjusted R-squared: 0.9262
## F-statistic: 51.19 on 1 and 3 DF, p-value: 0.005622The R packages and data we will be using.
The “dune” dataset we will use is included in the vegan package (an extensive collection of tools for multivariate analysis), but we will also be using some slightly modified and extended versions of it, which we will load as needed. We will also load some other packages to extend the range of analyses available to us.
##Q2
Question 2: Simple exercise for illustrating effects of different distance metrics. Should not be based on ordination or classification.
Simple exercise for illustrating effects of different distance metrics. Should not be based on ordination or classification.
#TBC, Martyn##Q3
Question 3: In this exercise we will use Principal Components Analysis on the dune meadow explanatory data. We will also explore the effects of different kinds of scaling on the ordination plot. NOTE that the scaling issues applies to all ordination techniques!
Compare the plots on the PCA with and without scaling (PCA plot 4) and discuss the differences. Is there any way to see that one is on scaled data and one is not? What is the difference between PCA plot 4 and PCA plot 5? Which one is best to use? Discuss the difference between the plots in plot PCA plot 6 and 7!
library(tidyverse)
library(vegan)
library(corrplot)
data("dune")
data("dune.env")
dummy_management <- as.data.frame(model.matrix( ~ Management - 1, data=dune.env ))
#add these to the dataset
dune.env.original <- dune.env #we keep a copy of the original version
dune.env <- dune.env %>% select(A1, Moisture, Manure, Use) %>% cbind(.,dummy_management)
dune.env$Moisture <- as.numeric(as.character(dune.env$Moisture)) #make numeric
dune.env$Manure <- as.numeric(as.character(dune.env$Manure))
dune.env$Use <- as.numeric(dune.env$Use)
#make column names shorter
dune.env <- dune.env %>% rename(BF = ManagementBF, HF = ManagementHF,
NM = ManagementNM, SF = ManagementSF)
## PCA plot 1:
# Start by looking at pairwise correlations among the variables.
dune.env_cor<-cor(dune.env, method = "kendall")
corrplot(dune.env_cor, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)
## PCA plot 2:
# Then make a PCA, with scaling of the data
env.pca <- rda(dune.env, scale = TRUE) #vegan uses the same function for PCA and RDA, just depends on if it is constrained or not.
biplot(env.pca) #plot the results using the default plot scaling which is "species"
env.pca #summarise results## Call: rda(X = dune.env, scale = TRUE)
##
## Inertia Rank
## Total 8
## Unconstrained 8 7
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## 2.7348 1.9195 1.3056 1.0477 0.5106 0.4141 0.0677summary(eigenvals(env.pca)) #see variance explained## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Eigenvalue 2.7348 1.9195 1.3056 1.048 0.51059 0.41413 0.067668
## Proportion Explained 0.3418 0.2399 0.1632 0.131 0.06382 0.05177 0.008458
## Cumulative Proportion 0.3418 0.5818 0.7450 0.876 0.93977 0.99154 1.000000## PCA plot 3:
# Continue with a PCA on the same data, but without scaling of the data
env.pca2 <- rda(dune.env, scale = FALSE) #vegan uses the same function for PCA and RDA, just depends on if it is constrained or not.
biplot(env.pca2)#plot the results
env.pca2 #summarise results## Call: rda(X = dune.env, scale = FALSE)
##
## Inertia Rank
## Total 11.59
## Unconstrained 11.59 7
## Inertia is variance
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## 6.436 2.314 1.993 0.435 0.255 0.126 0.035summary(eigenvals(env.pca2)) #see variance explained## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Eigenvalue 6.4360 2.3137 1.9929 0.43530 0.25488 0.12593 0.034957
## Proportion Explained 0.5551 0.1996 0.1719 0.03755 0.02198 0.01086 0.003015
## Cumulative Proportion 0.5551 0.7547 0.9266 0.96414 0.98612 0.99698 1.000000## PCA plot 4:
# Compare the PCAs with and without scaling of data
# Default scaling of the plots, which is "species"
par(mfrow = c(1,2))
biplot(env.pca, main = "PCA with scaling of data")
biplot(env.pca2, main = "PCA without scaling of data")
par(mfrow = c(1,1))
## PCA plot 5:
# Plots from the same PCA, but with plot scaling focused on sites and on species
par(mfrow = c(1,2))
biplot(env.pca, scaling = "sites", main = "Plot scaling on sites")
biplot(env.pca, scaling = "species", main = "Plot scaling on species")
par(mfrow = c(1,1))
## PCA plot 6:
# Then, a plot from the PCA on un-scaled data,
# but with plot scaling to mimic a PCA on scaled data (correlation = TRUE)
# Second plot is the PCA on scaled data
# Both plots with plot scaling focused on species
par(mfrow = c(1,2))
biplot(env.pca2, scaling = "species", correlation = TRUE, main = "PCA on un-sclaed data,\nplot scaling on correlations")
biplot(env.pca, scaling = "species", correlation = FALSE, main = "PCA on scaled data")
par(mfrow = c(1,1))
## PCA plot 7:
# Finally a plot from the PCA on un-scaled data,
# but with plot scaling to mimic a PCA on scaled data
# Second plot is the PCA on scaled data
# Both plots with plot scaling focused on sites
par(mfrow = c(1,2))
biplot(env.pca2, scaling = "sites", correlation = TRUE, main = "PCA on un-scaled data,\nplot scaling on correlations")
biplot(env.pca, scaling = "sites", main = "PCA on scaled data")
par(mfrow = c(1,1))##Q4
Question 4: Do a CA-ordination on the Dune Meadow species dataset.What are the results telling you?
Give a conceptual description on why objects/samples and descriptors/species to the left differ from objects and descriptors to the right, and those at the bottom from those at the top! (Plant ecologists may give a more detailed description, using their knowledge about the species in the dataset).
dune.ca <- cca(dune)
plot(dune.ca)
dune.ca## Call: cca(X = dune)
##
## Inertia Rank
## Total 2.115
## Unconstrained 2.115 19
## Inertia is scaled Chi-square
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.5360 0.4001 0.2598 0.1760 0.1448 0.1079 0.0925 0.0809
## (Showing 8 of 19 unconstrained eigenvalues)summary(eigenvals(dune.ca)) #proportion variance explained## Importance of components:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7
## Eigenvalue 0.5360 0.4001 0.2598 0.17598 0.14476 0.10791 0.09247
## Proportion Explained 0.2534 0.1892 0.1228 0.08319 0.06844 0.05102 0.04372
## Cumulative Proportion 0.2534 0.4426 0.5654 0.64858 0.71702 0.76804 0.81175
## CA8 CA9 CA10 CA11 CA12 CA13 CA14
## Eigenvalue 0.08091 0.07332 0.05630 0.04826 0.04125 0.03523 0.020529
## Proportion Explained 0.03825 0.03466 0.02661 0.02282 0.01950 0.01665 0.009705
## Cumulative Proportion 0.85000 0.88467 0.91128 0.93410 0.95360 0.97025 0.979955
## CA15 CA16 CA17 CA18 CA19
## Eigenvalue 0.014911 0.009074 0.007938 0.007002 0.003477
## Proportion Explained 0.007049 0.004290 0.003753 0.003310 0.001644
## Cumulative Proportion 0.987004 0.991293 0.995046 0.998356 1.000000##Q5
Question 5: Repeat exercise 4, but with DCA ordination instead.
Look at the eigenvalues, the length of gradient, the total variation and the ordination diagram. Explain the differences between results from CA and DCA.
dune.dca <- decorana(dune)
dune.dca ##
## Call:
## decorana(veg = dune)
##
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
##
## DCA1 DCA2 DCA3 DCA4
## Eigenvalues 0.5117 0.3036 0.12125 0.14267
## Decorana values 0.5360 0.2869 0.08136 0.04814
## Axis lengths 3.7004 3.1166 1.30055 1.47888#Detrended correspondence analysis (function decorana).
#Note that we do not get "variation explained" in the R implementation of DCA (and some other functions).
#Here, the developer explains why, "The total amount of variation is undefined in detrended
#correspondence analysis and therefore proportions from total are unknown and undefined.
#DCA is not a method for decomposition of variation, and therefore
#these proportions would not make sense either.
plot(dune.dca)
##Q6
Question 6: Add one new species with cover degree 9 at the site with the lowest number of species. Do a new DCA-ordination.
What happens with the eigenvalues and ordination diagram?
# calculate the number of species present in each plot
rowSums(dune != 0)## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 5 10 10 13 14 11 13 12 13 12 9 9 10 7 8 8 7 9 9 8# find the plot with the lowest number of species
which.min(rowSums(dune != 0))## 1
## 1# create a new data frame with a new species with cover degree 9 at site 1
NewSp <- c(9, rep(0,19))
SpeciesNewSp <- data.frame(dune, NewSp)
# run a new DCA with the extra species
DCANewSp <- decorana(SpeciesNewSp)
DCANewSp##
## Call:
## decorana(veg = SpeciesNewSp)
##
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
##
## DCA1 DCA2 DCA3 DCA4
## Eigenvalues 0.5229 0.3352 0.13870 0.09294
## Decorana values 0.5436 0.2636 0.08678 0.04037
## Axis lengths 3.4890 3.6245 1.49207 1.24554# plot the DCA
ordiplot(DCANewSp, type="text", main="DCA with new sp.")
You may be surprised that you haven’t got any total inertia values when printing decorana results, although in other software (e.g. CANOCO) these are available, together with percentage variance explained by particular axes. The reason for this is that DCA does not support the concept of total inertia values (also, it produces only four axes, i.e. four eigenvalues). See the comment in the code for the question 5.
##Q7
### Question 7: Effects from increased variation. Add one new sample to the data, then add a cover of 9 for the existing species Chenopodium album “Che_alb” in the new sample. Make a new DCA-ordination.
Next, remove the Che alb = 9, and add instead a cover of 5 for the exisiting species Leo_aut and make a new DCA.
What happens with the eigenvalues and ordination after adding Che alb and after adding Leo aut? Why did we choose to add these two species in the new sample?
# figure out the column number for Chenopodium album (abbreviated "Chenalbu")
which(colnames(dune)=="Che_alb")## integer(0)# create a new data frame with a new sample containing Chenopodium album with a cover of 9
NewSample <- c(rep(0,7), 9, rep(0,22))
SpeciesNewSample <- data.frame(rbind(dune, NewSample))
# run a new DCA with the extra sample
DCANewSample <- decorana(SpeciesNewSample)
DCANewSample##
## Call:
## decorana(veg = SpeciesNewSample)
##
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
##
## DCA1 DCA2 DCA3 DCA4
## Eigenvalues 0.9012 0.5109 0.3028 0.13538
## Decorana values 0.9020 0.5345 0.2856 0.09427
## Axis lengths 1.7784 3.6794 3.1985 1.51095# plot the DCA
ordiplot(DCANewSample, type="text", main="DCA with new sample Chenopodium album")
# figure out the column number for Leontodon autumnalis (abbreviated "Leo aut")
which(colnames(dune)=="Leo_aut")## integer(0)# create a new data frame with a new sample containing Leontodon autumnalis with a cover of 5
NewSample <- c(rep(0,15), 5, rep(0,10))
SpeciesNewSample <- data.frame(rbind(dune, NewSample))
# run a new DCA with the extra sample
DCANewSample <- decorana(SpeciesNewSample)
DCANewSample##
## Call:
## decorana(veg = SpeciesNewSample)
##
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
##
## DCA1 DCA2 DCA3 DCA4
## Eigenvalues 0.5125 0.3040 0.3030 0.12779
## Decorana values 0.5361 0.2869 0.2187 0.07207
## Axis lengths 3.7039 3.0498 0.5132 1.44059# plot the DCA
ordiplot(DCANewSample, type="text", main="DCA with new sample Leontodon autumnalis")
##Q8
Question 8: Non Metric Multidimensional Scaling
Compare the NMDS ordination with the ordinations obtained using CA and DCA! What is the stress value for the NMDS?
# NMDS on original data set with 2 reduced dimensions
NMDS <- metaMDS(dune, k=2)## Run 0 stress 0.1192678
## Run 1 stress 0.1812936
## Run 2 stress 0.1192678
## ... New best solution
## ... Procrustes: rmse 3.554456e-05 max resid 0.0001037297
## ... Similar to previous best
## Run 3 stress 0.2341478
## Run 4 stress 0.1808914
## Run 5 stress 0.180892
## Run 6 stress 0.1192678
## ... New best solution
## ... Procrustes: rmse 2.856648e-05 max resid 8.100342e-05
## ... Similar to previous best
## Run 7 stress 0.1192678
## ... Procrustes: rmse 4.60395e-05 max resid 0.0001366484
## ... Similar to previous best
## Run 8 stress 0.1808911
## Run 9 stress 0.1183186
## ... New best solution
## ... Procrustes: rmse 0.02027803 max resid 0.06500843
## Run 10 stress 0.1886532
## Run 11 stress 0.119269
## Run 12 stress 0.1183186
## ... New best solution
## ... Procrustes: rmse 2.017647e-05 max resid 6.602145e-05
## ... Similar to previous best
## Run 13 stress 0.1183186
## ... New best solution
## ... Procrustes: rmse 7.419144e-06 max resid 2.567168e-05
## ... Similar to previous best
## Run 14 stress 0.119268
## Run 15 stress 0.1183186
## ... Procrustes: rmse 2.117589e-05 max resid 5.801426e-05
## ... Similar to previous best
## Run 16 stress 0.1192678
## Run 17 stress 0.1183186
## ... Procrustes: rmse 2.214182e-05 max resid 6.882882e-05
## ... Similar to previous best
## Run 18 stress 0.119268
## Run 19 stress 0.119268
## Run 20 stress 0.1183186
## ... New best solution
## ... Procrustes: rmse 2.785698e-06 max resid 8.22521e-06
## ... Similar to previous best
## *** Solution reachedNMDS##
## Call:
## metaMDS(comm = dune, k = 2)
##
## global Multidimensional Scaling using monoMDS
##
## Data: dune
## Distance: bray
##
## Dimensions: 2
## Stress: 0.1183186
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'dune'# plot the NMDS
ordiplot(NMDS, type="text", main="NMDS")
#plot the CA
ordiplot(dune.ca, type="text", main="DCA")
#plot the DCA
ordiplot(dune.dca, type="text", main="DCA")
##Q9
Question 9: PCoA
data quality
#To follow, Martyn##Q10
Question 10: Data quality assessment
data quality
# To follow, Ulf##Q11
Question 11: Do an RDA and a CCA on the Dune Meadow data set, using the full set of environmental variables.
Compare the ordination diagrams and numerical results! What does the significance test tell you? When only looking at the results, which method explains the variation best? Taking the nature of the data into account, which method should be used?
# RDA
dune.rda <- rda(dune ~ . , dune.env.original)
summary(dune.rda)##
## Call:
## rda(formula = dune ~ A1 + Moisture + Management + Use + Manure, data = dune.env.original)
##
## Partitioning of variance:
## Inertia Proportion
## Total 84.12 1.0000
## Constrained 63.21 0.7513
## Unconstrained 20.92 0.2487
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7
## Eigenvalue 22.3955 16.2076 7.03891 4.0380 3.7602 2.60874 2.16693
## Proportion Explained 0.2662 0.1927 0.08367 0.0480 0.0447 0.03101 0.02576
## Cumulative Proportion 0.2662 0.4589 0.54256 0.5906 0.6353 0.66627 0.69203
## RDA8 RDA9 RDA10 RDA11 RDA12 PC1 PC2
## Eigenvalue 1.80327 1.40421 0.91739 0.581545 0.283927 6.62687 4.30914
## Proportion Explained 0.02144 0.01669 0.01091 0.006913 0.003375 0.07878 0.05122
## Cumulative Proportion 0.71346 0.73015 0.74106 0.747973 0.751348 0.83012 0.88135
## PC3 PC4 PC5 PC6 PC7
## Eigenvalue 3.54913 2.54648 2.34027 0.9335 0.612096
## Proportion Explained 0.04219 0.03027 0.02782 0.0111 0.007276
## Cumulative Proportion 0.92354 0.95381 0.98163 0.9927 1.000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7
## Eigenvalue 22.3955 16.2076 7.0389 4.03795 3.76017 2.60874 2.16693
## Proportion Explained 0.3543 0.2564 0.1114 0.06389 0.05949 0.04127 0.03428
## Cumulative Proportion 0.3543 0.6107 0.7221 0.78600 0.84549 0.88676 0.92105
## RDA8 RDA9 RDA10 RDA11 RDA12
## Eigenvalue 1.80327 1.40421 0.91739 0.581545 0.283927
## Proportion Explained 0.02853 0.02222 0.01451 0.009201 0.004492
## Cumulative Proportion 0.94958 0.97179 0.98631 0.995508 1.000000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 6.322924
##
##
## Species scores
##
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## Achimill 0.58181 0.08704 -0.105553 -0.043694 -0.422290 -0.062221
## Agrostol -1.16975 -0.94792 0.060309 0.334142 0.030721 0.222829
## Airaprae -0.11716 0.20786 0.018154 -0.180185 -0.050376 0.041993
## Alopgeni -0.50957 -1.41677 -0.322262 -0.116416 -0.013318 -0.300469
## Anthodor 0.39549 0.45832 -0.411832 -0.136384 -0.336365 0.263272
## Bellpere 0.36549 -0.10133 0.162631 -0.160330 -0.370665 0.137716
## Bromhord 0.53405 -0.19418 0.161217 -0.048422 -0.585999 -0.030597
## Chenalbu -0.03146 -0.03881 -0.034281 0.004236 -0.017102 0.003806
## Cirsarve 0.02198 -0.09707 0.075983 0.009750 0.017309 0.069692
## Comapalu -0.19434 0.09464 0.045000 0.157251 -0.021772 -0.021112
## Eleopalu -1.07946 0.08439 0.006165 0.698662 -0.258796 0.127090
## Elymrepe 0.57194 -0.67735 0.299352 -0.407119 0.233601 0.408690
## Empenigr -0.07889 0.08081 0.032442 -0.076400 -0.005971 0.030773
## Hyporadi -0.12994 0.36832 0.135169 -0.148462 0.009242 -0.069277
## Juncarti -0.46506 -0.11361 -0.043146 0.023106 0.051743 0.032617
## Juncbufo -0.13148 -0.45491 -0.399997 -0.320717 0.288602 -0.369510
## Lolipere 1.57085 -0.39016 0.727211 0.554451 0.340886 -0.082511
## Planlanc 0.89986 0.57656 -0.564648 0.273286 0.038960 0.120162
## Poaprat 0.98798 -0.44183 0.418705 0.121187 0.153703 0.036591
## Poatriv 0.74253 -1.39782 -0.559135 0.008607 -0.428312 0.162203
## Ranuflam -0.55236 0.08108 0.004173 0.185963 -0.131421 0.083262
## Rumeacet 0.57938 0.06076 -0.940669 0.137799 0.309740 0.163586
## Sagiproc -0.21299 -0.45646 0.011857 -0.174011 0.208064 -0.216092
## Salirepe -0.28603 0.50726 0.148922 -0.456522 0.030091 0.164904
## Scorautu 0.33487 0.50344 0.086424 -0.148075 -0.258282 -0.288324
## Trifprat 0.42084 0.16346 -0.528815 0.260068 0.196986 0.120890
## Trifrepe 0.39232 0.08889 -0.146355 0.260759 -0.280784 -0.613465
## Vicilath 0.11379 0.13700 0.111709 0.049534 0.011107 -0.146626
## Bracruta -0.01603 0.41968 -0.395909 0.008729 0.367890 -0.142926
## Callcusp -0.52015 0.21794 0.061946 0.231494 -0.019057 0.025981
##
##
## Site scores (weighted sums of species scores)
##
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## 1 1.01596 -0.1799 2.59463 -0.11443 2.05964 2.1377
## 2 1.72853 -1.0870 1.31739 -0.89469 -3.45007 -1.0625
## 3 0.61315 -2.4671 1.08784 -0.06821 0.64265 0.3666
## 4 0.12048 -2.1589 1.71771 -0.07977 0.48648 1.8652
## 5 1.68284 0.4191 -2.55443 -0.80052 -1.22481 3.3283
## 6 2.08298 1.2151 -3.01397 2.45080 1.69576 -0.3949
## 7 1.89757 0.4084 -1.15310 1.15696 0.19469 0.1648
## 8 -0.79508 -1.2469 0.59088 1.85859 0.03913 -0.7791
## 9 0.08193 -1.7029 -0.41962 -2.00241 1.90642 0.0328
## 10 1.95184 0.6210 0.49322 0.87955 -3.21475 -1.7550
## 11 1.07928 1.3401 1.77536 0.93077 2.20291 -2.9068
## 12 -1.19429 -1.6563 -2.36651 -1.91407 1.14492 -3.7773
## 13 -0.82017 -2.3013 -1.46058 -1.16242 -1.40595 -0.6812
## 14 -1.69462 1.0179 0.55226 2.12686 -1.38338 -1.1655
## 15 -1.92197 1.0316 0.18669 1.38332 0.21043 0.9451
## 16 -2.81781 -0.4853 -0.45260 2.98066 -0.24307 2.1409
## 17 -0.26124 1.7408 -0.04449 -2.00579 -0.85795 1.4187
## 18 0.25044 1.9245 0.47979 -1.27429 1.11822 -1.1395
## 19 -0.64055 2.2851 0.18205 -3.72278 -0.40927 -1.1779
## 20 -2.35927 1.2820 0.48747 0.27188 0.48801 2.4397
##
##
## Site constraints (linear combinations of constraining variables)
##
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## 1 0.87487 -0.90608 2.36103 -0.66159 2.53531 1.55448
## 2 1.91194 -0.62145 1.07037 -0.88335 -3.00393 -2.04867
## 3 0.43540 -1.95933 1.52026 0.20428 0.30748 1.42187
## 4 0.43933 -1.94048 1.51888 0.19490 0.34600 1.39312
## 5 1.65691 0.51332 -2.47861 0.06329 -1.37887 3.19410
## 6 2.10890 1.12091 -3.08978 1.58700 1.84982 -0.26074
## 7 1.48330 -0.04811 -0.36778 1.16790 0.69199 -0.12568
## 8 -0.56422 -1.25593 0.05259 1.83632 -0.90431 0.49765
## 9 0.26534 -1.23731 -0.66664 -1.99108 2.35256 -0.95341
## 10 1.52623 -0.09659 0.82501 -0.36238 -2.76909 -0.83620
## 11 1.32148 1.59213 1.69059 2.16136 1.31111 -2.83941
## 12 -1.37770 -2.12186 -2.11948 -1.92541 0.69878 -2.79106
## 13 -1.25782 -1.55141 -1.37052 0.16936 -0.68373 0.15215
## 14 -2.37447 1.51577 0.72709 2.22806 0.23064 -0.39953
## 15 -1.51025 0.37596 0.17244 0.91534 -0.66586 -0.02248
## 16 -2.19675 -0.76968 -0.78968 1.66022 -0.51916 0.32129
## 17 0.02328 1.73218 -0.60987 -1.31102 -0.82795 -0.08329
## 18 0.38018 2.38966 0.25985 -1.97999 0.59091 0.65301
## 19 -1.57691 1.61530 0.64851 -1.52722 -0.11937 0.61515
## 20 -1.56904 1.65300 0.64575 -1.54597 -0.04232 0.55766
##
##
## Biplot scores for constraining variables
##
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## A1 -0.5361 0.06882 -0.218452 0.31811 -0.146216 0.05232
## Moisture.L -0.9078 -0.11664 -0.056742 0.05789 -0.136471 -0.01179
## Moisture.Q -0.1012 0.43397 -0.036128 0.40058 -0.008317 0.14269
## Moisture.C -0.1978 0.01609 0.364691 0.19843 -0.452618 0.35989
## ManagementHF 0.4043 -0.07409 -0.534962 0.21752 0.213258 0.19208
## ManagementNM -0.5114 0.71630 0.142286 -0.24855 -0.064357 0.10191
## ManagementSF -0.2379 -0.71375 0.086470 -0.02765 0.207180 0.15834
## Use.L -0.1398 -0.25485 0.123795 0.81530 0.138325 -0.26210
## Use.Q -0.1173 0.54889 -0.004563 0.02884 -0.099898 0.04341
## Manure.L 0.2718 -0.78470 0.068253 0.26784 0.140391 0.28325
## Manure.Q -0.3913 0.30047 0.585951 -0.19204 0.183926 0.39222
## Manure.C 0.5204 -0.16720 0.322839 -0.18763 0.228598 -0.20896
##
##
## Centroids for factor constraints
##
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## Moisture1 1.3911 0.57720 -0.079190 0.20780 0.3709 0.01816
## Moisture2 0.6061 -0.56605 0.813569 -0.31856 -0.7359 0.47388
## Moisture4 -0.5562 -1.67959 -1.393063 -1.95824 1.5257 -1.87223
## Moisture5 -1.5785 0.22614 0.012311 0.53373 -0.3863 0.24598
## ManagementBF 1.5865 0.29136 1.195323 0.30521 -1.4873 -1.90809
## ManagementHF 0.9900 -0.18143 -1.310046 0.53268 0.5222 0.47038
## ManagementNM -1.1045 1.54698 0.307294 -0.53680 -0.1390 0.22008
## ManagementSF -0.5138 -1.54147 0.186749 -0.05971 0.4474 0.34198
## UseHayfield 0.1009 0.93851 -0.196574 -1.23634 -0.3134 0.44957
## UseHaypastu 0.2031 -0.95047 0.007901 -0.04993 0.1730 -0.07517
## UsePasture -0.4661 0.20684 0.262562 1.81077 0.1621 -0.50914
## Manure1 1.0377 0.08608 0.616318 -0.06403 0.2982 -1.54300
## Manure2 1.0750 -0.27727 -1.654376 -0.28962 -0.4585 -0.47659
## Manure3 -0.6339 -0.90628 -0.618849 1.20845 -0.3538 0.21135
## Manure4 0.5832 -1.60196 1.800060 -0.08747 1.0629 1.45649anova(dune.rda)# CCA
dune.cca <- cca(dune ~ . , dune.env.original)
summary(dune.cca)##
## Call:
## cca(formula = dune ~ A1 + Moisture + Management + Use + Manure, data = dune.env.original)
##
## Partitioning of scaled Chi-square:
## Inertia Proportion
## Total 2.1153 1.0000
## Constrained 1.5032 0.7106
## Unconstrained 0.6121 0.2894
##
## Eigenvalues, and their contribution to the scaled Chi-square
##
## Importance of components:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7
## Eigenvalue 0.4671 0.3410 0.17606 0.15317 0.09528 0.07027 0.05887
## Proportion Explained 0.2208 0.1612 0.08323 0.07241 0.04504 0.03322 0.02783
## Cumulative Proportion 0.2208 0.3821 0.46529 0.53770 0.58275 0.61596 0.64379
## CCA8 CCA9 CCA10 CCA11 CCA12 CA1 CA2
## Eigenvalue 0.04993 0.03183 0.02596 0.02282 0.010818 0.2724 0.10876
## Proportion Explained 0.02360 0.01505 0.01227 0.01079 0.005114 0.1288 0.05142
## Cumulative Proportion 0.66740 0.68245 0.69472 0.70551 0.710627 0.8394 0.89081
## CA3 CA4 CA5 CA6 CA7
## Eigenvalue 0.08975 0.06305 0.03489 0.02529 0.017984
## Proportion Explained 0.04243 0.02981 0.01649 0.01196 0.008502
## Cumulative Proportion 0.93324 0.96305 0.97954 0.99150 1.000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7
## Eigenvalue 0.4671 0.3410 0.1761 0.1532 0.09528 0.07027 0.05887
## Proportion Explained 0.3108 0.2269 0.1171 0.1019 0.06339 0.04674 0.03916
## Cumulative Proportion 0.3108 0.5376 0.6548 0.7567 0.82005 0.86679 0.90595
## CCA8 CCA9 CCA10 CCA11 CCA12
## Eigenvalue 0.04993 0.03183 0.02596 0.02282 0.010818
## Proportion Explained 0.03322 0.02118 0.01727 0.01518 0.007197
## Cumulative Proportion 0.93917 0.96035 0.97762 0.99280 1.000000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
##
##
## Species scores
##
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6
## Achimill 0.857707 0.49484 -0.045074 -0.057950 0.670679 0.056946
## Agrostol -0.750987 -0.51041 -0.001514 -0.006522 0.110328 -0.113279
## Airaprae -0.786606 1.84744 0.777309 1.068417 0.558936 -0.142093
## Alopgeni -0.296764 -0.98554 0.015875 0.454024 0.295335 0.124614
## Anthodor 0.366346 1.04925 -0.217724 0.412435 0.563903 -0.226744
## Bellpere 0.707119 0.06620 0.500239 -0.369709 0.592272 -0.330213
## Bromhord 0.853359 -0.02616 0.335626 -0.427024 0.865565 -0.093885
## Chenalbu -0.672293 -1.27568 -0.417041 0.994310 0.640906 0.521336
## Cirsarve 0.386092 -0.93251 0.963978 -0.468823 0.239146 -0.908823
## Comapalu -2.015388 0.10446 -1.185954 -1.972940 -0.083546 -0.219456
## Eleopalu -1.515611 -0.03069 -0.455387 -0.508862 0.111183 -0.132850
## Elymrepe 0.610453 -0.66152 0.628614 -0.124262 -0.491387 -0.510855
## Empenigr -1.378845 1.51426 1.387645 0.973499 -0.434196 -0.642379
## Hyporadi -0.655401 1.49855 0.738512 0.362909 -0.104564 0.310460
## Juncarti -0.948447 -0.19092 -0.028958 0.148431 -0.292973 -0.155408
## Juncbufo -0.064392 -1.01588 -0.297937 1.118684 -0.196979 0.823725
## Lolipere 0.607215 -0.09885 0.262228 -0.409402 -0.312573 0.082660
## Planlanc 0.800409 0.77629 -0.606615 -0.038048 -0.021902 0.152465
## Poaprat 0.471703 -0.13797 0.252857 -0.205726 -0.176841 0.061944
## Poatriv 0.303929 -0.50098 -0.033488 0.157680 0.180752 -0.110423
## Ranuflam -1.384091 0.11712 -0.129524 -0.120128 0.041536 -0.207804
## Rumeacet 0.791036 0.13052 -1.100806 0.537567 -0.416656 -0.396323
## Sagiproc -0.290759 -0.52139 0.288191 0.373446 -0.109281 0.193977
## Salirepe -0.889942 1.59655 1.182347 0.817665 -0.457263 -0.271421
## Scorautu -0.009838 0.37963 0.098333 -0.075560 0.001781 0.125526
## Trifprat 1.054359 0.56735 -1.455764 0.398295 -0.539708 -0.564846
## Trifrepe 0.014698 0.11386 -0.184205 -0.244273 0.089510 0.207091
## Vicilath 0.511432 0.88591 0.273760 -1.050861 -0.335122 1.591159
## Bracruta -0.147674 0.27946 -0.137684 0.084272 -0.210685 0.147748
## Callcusp -1.773372 0.55553 -0.188876 -0.426816 -0.097141 -0.008786
##
##
## Site scores (weighted averages of species scores)
##
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6
## 1 1.19460 -0.71633 1.656429 -1.4249 -2.23242 -1.09183
## 2 0.86805 -0.35777 0.934868 -0.9242 1.33133 -0.48738
## 3 0.33794 -1.04078 0.824936 -0.1786 -0.04887 -0.42012
## 4 0.25095 -0.99294 1.237259 -0.4411 0.35593 -1.14729
## 5 1.11991 0.45932 -1.022763 0.2895 0.27912 -1.76102
## 6 0.99305 0.73388 -2.001441 0.3758 -1.04277 -0.73162
## 7 1.03098 0.34363 -1.083231 0.2079 -0.14394 0.29945
## 8 -0.66671 -0.71037 0.004385 -0.1742 -0.10923 0.07552
## 9 0.09269 -1.09341 0.218954 0.9201 -1.32966 0.09022
## 10 0.95315 0.58996 0.146744 -0.9667 1.81728 0.84400
## 11 0.47318 0.74856 0.535395 -1.0620 -1.51667 3.25324
## 12 -0.27934 -1.30695 -0.512852 2.0022 0.26928 1.96951
## 13 -0.37400 -1.45815 -0.153034 1.4764 1.08084 1.21481
## 14 -2.04173 0.23744 -1.448971 -2.6038 0.41722 -0.22568
## 15 -1.93799 -0.04255 -1.367867 -1.8081 -0.32644 -0.85946
## 16 -1.91272 -0.56130 -0.881308 -0.5388 0.40537 -0.81739
## 17 0.33533 2.74717 0.494859 1.7252 3.00069 0.15330
## 18 0.26805 1.23029 0.753438 -0.3018 -1.09614 1.50399
## 19 -0.75573 2.53794 2.225936 2.4221 0.13516 0.14348
## 20 -2.03905 0.80938 0.463176 0.2188 -1.15999 -1.41953
##
##
## Site constraints (linear combinations of constraining variables)
##
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6
## 1 0.7577 -1.0533 2.075818 -0.62357 -2.6193 -0.57163
## 2 1.0307 -0.2511 0.942759 -0.72795 1.0939 0.10633
## 3 0.3825 -0.9571 0.943652 -0.50806 0.2566 -0.92249
## 4 0.3861 -0.9325 0.963978 -0.46882 0.2391 -0.90882
## 5 1.0281 0.3743 -1.180481 -0.15376 0.1863 -2.04984
## 6 1.0753 0.8100 -1.860152 0.77287 -0.9596 -0.47288
## 7 1.0282 0.1537 -0.720078 0.01393 -0.2159 0.69023
## 8 -0.8347 -0.6324 -0.367053 -0.18633 0.2121 -0.93866
## 9 0.2553 -0.9868 0.226845 1.11630 -1.5671 0.68393
## 10 0.8517 0.5884 0.543136 -0.45186 1.9874 0.90302
## 11 0.3961 0.6106 -0.007613 -2.01145 -1.4337 2.39468
## 12 -0.4745 -1.4349 -0.522321 1.76672 0.5542 1.25706
## 13 -0.6723 -1.2757 -0.417041 0.99431 0.6409 0.52134
## 14 -2.3490 0.6722 -1.087969 -2.00082 -0.2983 -0.03446
## 15 -1.6817 -0.4633 -1.283939 -1.94506 0.1312 -0.40445
## 16 -1.4075 -0.6080 -0.607258 0.19306 0.5431 0.63171
## 17 0.1018 2.3472 -0.138197 1.21079 2.0486 0.60834
## 18 0.4019 1.7339 0.567131 0.27132 -0.4606 0.67225
## 19 -1.3788 1.5143 1.387645 0.97350 -0.4342 -0.64238
## 20 -1.3717 1.5635 1.428297 1.05197 -0.4691 -0.61505
##
##
## Biplot scores for constraining variables
##
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6
## A1 -0.5534 -0.14431 -0.51439 -0.29750 0.10556 -0.18207
## Moisture.L -0.9117 -0.18833 0.03981 0.20831 0.10483 -0.03496
## Moisture.Q -0.2030 0.34831 -0.29281 -0.20645 -0.28387 -0.13407
## Moisture.C -0.2100 0.10223 0.30648 -0.29962 0.52384 -0.28198
## ManagementHF 0.3604 -0.01771 -0.54879 0.22226 -0.32807 -0.29032
## ManagementNM -0.5914 0.64415 0.15809 -0.01577 -0.04877 -0.07689
## ManagementSF -0.1243 -0.67548 0.21199 0.15443 0.10327 -0.01928
## Use.L -0.1943 -0.34773 -0.33160 -0.42673 -0.10227 0.22872
## Use.Q -0.1889 0.51516 -0.02141 -0.01123 -0.05134 0.16126
## Manure.L 0.3295 -0.69510 0.03992 -0.02906 0.04350 -0.26625
## Manure.Q -0.4229 0.22292 0.52833 -0.20228 -0.12773 -0.23201
## Manure.C 0.4853 -0.11991 0.31222 -0.19850 -0.15170 0.17650
##
##
## Centroids for factor constraints
##
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6
## Moisture1 0.86961 0.39190 -0.30729 -0.271184 -0.42482 0.0229
## Moisture2 0.49529 -0.13803 0.71613 -0.298514 0.95955 -0.2087
## Moisture4 -0.07641 -1.19047 -0.11369 1.411949 -0.60287 0.9444
## Moisture5 -1.31323 0.06247 -0.07829 0.008202 0.07168 -0.2261
## ManagementBF 0.79133 0.29313 0.53596 -0.977522 0.73099 1.0250
## ManagementHF 0.53656 -0.02637 -0.81693 0.330865 -0.48837 -0.4322
## ManagementNM -1.11223 1.21134 0.29730 -0.029662 -0.09172 -0.1446
## ManagementSF -0.19081 -1.03722 0.32552 0.237132 0.15857 -0.0296
## UseHayfield 0.08045 0.76460 0.35627 0.470345 0.07734 -0.1392
## UseHaypastu 0.22445 -0.61214 0.02544 0.013349 0.06101 -0.1916
## UsePasture -0.48763 -0.02094 -0.53183 -0.668114 -0.20870 0.5131
## Manure1 0.51300 0.02905 0.27896 -0.315483 -0.22424 1.2323
## Manure2 0.72919 -0.03448 -0.70675 0.367543 0.16244 -0.3713
## Manure3 -0.41706 -0.55693 -0.52936 0.221143 0.26656 0.1926
## Manure4 0.44964 -0.96318 1.15039 -0.511103 -0.25361 -0.8552anova(dune.cca)#plots
ordiplot(dune.rda, type="text", main="RDA")
ordiplot(dune.cca, type="text", main="CCA")
##Q12
Question 12: Use the results from the CCA ordination in exercise 13 to look at and interpret the variance inflation factor, and to look for warnings about outliers in the explanatory data!
Look at and interpret the inflation factors for each environmental variable! Are there any extreme values?
# VIF for CCA
dune.cca <- cca(dune ~ . , dune.env)
vif.cca <- vif.cca(dune.cca)
vif.cca## A1 Moisture Manure Use BF HF NM SF
## 1.523245 1.733303 7.878266 1.683073 4.685483 2.945088 7.490125 NAFunction vif.cca gives the variance inflation factors for each constraint or contrast in factor constraints. In partial ordination, conditioning variables are analysed together with constraints. Variance inflation is a diagnostic tool to identify useless constraints. A common rule is that values over 10 indicate redundant constraints.
##Q13
Question 13: Do a CA on the Dune Meadow data set, with the full set of environmental variables included in the analysis, and compare with a CCA on the same data sets.
Compare the constrained and unconstrained ordinations! What is the difference between the methods? When could an unconstrained ordination with explanatory variables be a good option?
# CA with env. variables
fit <- envfit(dune.ca, dune.env, perm=999)
fit##
## ***VECTORS
##
## CA1 CA2 r2 Pr(>r)
## A1 0.99816 0.06061 0.3104 0.046 *
## Moisture 0.95470 0.29756 0.6802 0.001 ***
## Manure -0.16779 -0.98582 0.4848 0.003 **
## Use 0.73455 -0.67855 0.1616 0.175
## BF -0.98156 -0.19113 0.1122 0.311
## HF -0.79401 -0.60790 0.1052 0.366
## NM 0.41219 0.91110 0.7054 0.001 ***
## SF 0.44670 -0.89468 0.2428 0.066 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999ordiplot(dune.ca, type="text", main="CA with env. variables fitted after")
plot(fit)
#CCA
dune.cca <- cca(dune ~ . , dune.env.original)
dune.cca## Call: cca(formula = dune ~ A1 + Moisture + Management + Use + Manure,
## data = dune.env.original)
##
## Inertia Proportion Rank
## Total 2.1153 1.0000
## Constrained 1.5032 0.7106 12
## Unconstrained 0.6121 0.2894 7
## Inertia is scaled Chi-square
## Some constraints were aliased because they were collinear (redundant)
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10 CCA11
## 0.4671 0.3410 0.1761 0.1532 0.0953 0.0703 0.0589 0.0499 0.0318 0.0260 0.0228
## CCA12
## 0.0108
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7
## 0.27237 0.10876 0.08975 0.06305 0.03489 0.02529 0.01798ordiplot(dune.cca, type="text", main="CA")
#Q14 ### Question 14: Site scores can be calculated in two ways; Weighted Average and Linear Combination. Here you will do both and compare the results. Both options are available in most software, but the default varies.
What is the main difference between the two ways. How to know which method to select? What does the display = “cn” in the last plots mean?
#### Illustration of difference between WA and LC scores
data("dune")
data("dune.env")
dune.env$Moisture<-as.numeric(dune.env$Moisture) #Set Moisture to numeric for the illustration of WA and LC
# Make a CCA with Management and Moisture as explanatory variables
cca.dune<-cca(dune~Management + Moisture, dune.env)
# Look at the result
anova(cca.dune)summary(eigenvals(cca.dune))## Importance of components:
## CCA1 CCA2 CCA3 CCA4 CA1 CA2 CA3
## Eigenvalue 0.4355 0.2722 0.10304 0.04575 0.3709 0.16307 0.13787
## Proportion Explained 0.2059 0.1287 0.04871 0.02163 0.1754 0.07709 0.06518
## Cumulative Proportion 0.2059 0.3345 0.38326 0.40488 0.5803 0.65734 0.72252
## CA4 CA5 CA6 CA7 CA8 CA9 CA10
## Eigenvalue 0.12969 0.09557 0.07909 0.06565 0.05511 0.04546 0.03917
## Proportion Explained 0.06131 0.04518 0.03739 0.03104 0.02605 0.02149 0.01852
## Cumulative Proportion 0.78383 0.82901 0.86640 0.89744 0.92349 0.94498 0.96350
## CA11 CA12 CA13 CA14 CA15
## Eigenvalue 0.02712 0.020921 0.016962 0.008008 0.004208
## Proportion Explained 0.01282 0.009891 0.008019 0.003786 0.001989
## Cumulative Proportion 0.97632 0.986206 0.994225 0.998011 1.000000# Make two plots, one with weighted average (WA) and one with Linear Combination (LC) to calculate the site scores.
par(mfrow= c(1,2)) # To get two plots beside each other
plot(cca.dune, type = "n", display = c("wa"), scaling = "sites", main = "WA", pch = 3)
points(cca.dune, display = c("wa"), scaling = "sites", pch = 3)
text(cca.dune, display = "wa", scaling = "sites", pos = 2)
ordispider(cca.dune, dune.env$Management, col = "black", label= TRUE, scaling = "sites", display = "wa")
plot(cca.dune, type = "n", display = c("lc"), scaling = "sites", main = "LC", pch = 3)
points(cca.dune, display = c("lc"), scaling = "sites", col = "red", pch = 19)
text(cca.dune, display = "lc", col = "red", scaling = "sites", pos = 3)
ordispider(cca.dune, dune.env$Management, col = "red", label= TRUE, scaling = "sites", display = "lc")
par(mfrow= c(1,1))
# Same as above, but here with sites and species
par(mfrow= c(1,2))
plot(cca.dune, display = c("sp", "wa", "cn"), scaling = "species", main = "Weighted average")
plot(cca.dune, display = c("sp", "lc", "cn"), scaling = "species", main = "Linear combination")
par(mfrow= c(1,1))
##Q15
Question 15: Next, you should investigate whether the four different management regimes (SF standard farming, BF biodynamical farming, HF hobby farming, and NM nature management) have any effect on the species distribution. This is done by using these four management types as environmental variables in a constrained ordination.
Inspect and interpret the results of the permutations and graph! Compare with the CCA on the full set of explanatory data in exercise 13.
To investigate if the vegetational pattern is an effect of management, or if the pattern is random, Monte Carlo permutations will be used.
# CCA using only Management as the environmental variable
dune.ccaMN <- cca(dune ~ Management , dune.env.original)
summary(dune.ccaMN)##
## Call:
## cca(formula = dune ~ Management, data = dune.env.original)
##
## Partitioning of scaled Chi-square:
## Inertia Proportion
## Total 2.1153 1.0000
## Constrained 0.6038 0.2855
## Unconstrained 1.5114 0.7145
##
## Eigenvalues, and their contribution to the scaled Chi-square
##
## Importance of components:
## CCA1 CCA2 CCA3 CA1 CA2 CA3 CA4
## Eigenvalue 0.3186 0.18247 0.10274 0.4474 0.20300 0.16301 0.13457
## Proportion Explained 0.1506 0.08626 0.04857 0.2115 0.09597 0.07706 0.06362
## Cumulative Proportion 0.1506 0.23690 0.28547 0.4970 0.59294 0.67000 0.73362
## CA5 CA6 CA7 CA8 CA9 CA10 CA11
## Eigenvalue 0.12940 0.09494 0.07904 0.06526 0.05004 0.04321 0.0387
## Proportion Explained 0.06117 0.04488 0.03737 0.03085 0.02366 0.02043 0.0183
## Cumulative Proportion 0.79479 0.83968 0.87704 0.90790 0.93155 0.95198 0.9703
## CA12 CA13 CA14 CA15 CA16
## Eigenvalue 0.02385 0.017729 0.009172 0.007959 0.004157
## Proportion Explained 0.01128 0.008382 0.004336 0.003763 0.001965
## Cumulative Proportion 0.98155 0.989936 0.994272 0.998035 1.000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## CCA1 CCA2 CCA3
## Eigenvalue 0.3186 0.1825 0.1027
## Proportion Explained 0.5277 0.3022 0.1701
## Cumulative Proportion 0.5277 0.8299 1.0000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
##
##
## Species scores
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## Achimill 0.19825 0.72622 0.300280 0.536978 -0.435913 0.64609
## Agrostol -0.10058 -0.73046 0.007614 -0.681142 0.217728 0.03744
## Airaprae -1.87848 -0.05503 -0.068523 2.450771 1.553306 0.96472
## Alopgeni 0.55258 -0.78563 0.076123 -0.233871 0.552230 -0.21144
## Anthodor -0.39419 0.50501 -0.225911 1.160620 0.171184 0.93294
## Bellpere 0.13512 0.17366 0.564850 0.257160 -0.708045 -0.32204
## Bromhord 0.49092 0.60069 0.653072 0.176934 -0.436949 0.13777
## Chenalbu 0.56015 -1.38599 0.350857 0.002875 0.943749 -0.06831
## Cirsarve 0.56015 -1.38599 0.350857 0.405234 -0.657911 -0.47130
## Comapalu -1.87848 -0.05503 -0.068523 -1.746112 -0.512810 0.93566
## Eleopalu -0.70823 -0.37016 -0.119191 -1.553248 0.282390 0.66025
## Elymrepe 0.53962 -0.19035 -0.046524 0.198907 -0.499457 -0.61438
## Empenigr -1.87848 -0.05503 -0.068523 2.309910 2.212605 -0.08448
## Hyporadi -1.36521 0.25217 0.325575 1.904376 1.495023 0.41829
## Juncarti -0.38902 0.03085 -0.512153 -1.205380 0.574275 -0.35684
## Juncbufo 0.55931 -0.45216 -0.375981 0.003634 0.784787 -0.51755
## Lolipere 0.43543 0.19968 0.213581 0.216193 -0.477291 -0.18621
## Planlanc 0.06038 0.66342 -0.325865 0.556636 -0.525653 0.42479
## Poaprat 0.32406 0.11985 0.096677 0.195454 -0.309623 -0.29579
## Poatriv 0.53695 -0.14144 -0.012649 0.046730 -0.006836 -0.07777
## Ranuflam -0.83362 -0.33640 -0.113762 -1.301757 0.264387 0.17231
## Rumeacet 0.55852 0.41250 -1.048979 0.325485 -0.319841 0.40174
## Sagiproc 0.18110 -0.51035 0.108393 0.359930 0.783570 -0.45162
## Salirepe -1.87848 -0.05503 -0.068523 0.123760 0.058476 -0.95975
## Scorautu -0.32939 0.26009 0.120989 0.263833 0.097401 -0.12059
## Trifprat 0.55832 0.63731 -1.223958 0.430250 -0.704454 0.82922
## Trifrepe -0.04953 0.33644 0.186947 -0.090075 -0.056136 0.12455
## Vicilath -0.14617 0.98179 1.261557 0.241368 -0.211654 -0.51063
## Bracruta -0.30221 -0.01389 -0.078786 -0.087945 -0.022522 -0.11640
## Callcusp -1.14689 -0.45432 0.057291 -1.645373 -0.097690 0.92131
##
##
## Site scores (weighted averages of species scores)
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## 1 1.35561 0.4747 1.06564 0.823027 -2.26986 -0.78064
## 2 0.86419 0.8142 1.87245 -0.198170 -0.21094 -0.45453
## 3 0.94050 -0.9209 0.88921 0.342743 -0.72769 -0.44672
## 4 0.77228 -1.1875 1.30265 0.405234 -0.65791 -0.47130
## 5 0.69269 1.5459 -1.41512 0.521056 -0.79020 0.65898
## 6 0.54544 1.7488 -2.38761 0.404007 -0.73194 1.05302
## 7 0.76144 1.5049 -1.05537 0.405051 -0.54999 0.43997
## 8 0.01261 -1.0636 -0.08646 -1.029299 1.29069 -0.66249
## 9 0.76176 -0.7696 -0.94225 -0.400662 0.94009 -1.66619
## 10 0.41252 1.8518 1.64522 0.094693 -0.26410 0.77515
## 11 -0.11176 1.2962 1.56524 0.132854 0.63174 -0.44504
## 12 0.73786 -1.7731 -0.58173 0.207790 1.17765 -0.18459
## 13 0.75284 -1.9690 0.09623 0.002875 0.94375 -0.06831
## 14 -1.85731 -1.0193 0.31713 -1.775732 -0.79413 1.57589
## 15 -1.69911 -1.1108 -0.99568 -1.716492 -0.23149 0.29543
## 16 -1.00521 -2.2417 -0.74574 -1.640220 0.82453 1.87405
## 17 -1.64906 2.1313 -0.06623 2.662063 0.56436 2.53853
## 18 -0.81907 1.1779 0.87653 0.605071 -1.84599 -1.92760
## 19 -2.61902 0.6190 0.42849 2.309910 2.21260 -0.08448
## 20 -2.32115 -1.3309 -1.00044 -1.476716 -0.09132 -0.90420
##
##
## Site constraints (linear combinations of constraining variables)
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## 1 0.5601 -1.38599 0.35086 0.823027 -2.26986 -0.78064
## 2 0.4313 1.32739 1.70492 -0.198170 -0.21094 -0.45453
## 3 0.5601 -1.38599 0.35086 0.342743 -0.72769 -0.44672
## 4 0.5601 -1.38599 0.35086 0.405234 -0.65791 -0.47130
## 5 0.5583 0.63731 -1.22396 0.521056 -0.79020 0.65898
## 6 0.5583 0.63731 -1.22396 0.404007 -0.73194 1.05302
## 7 0.5583 0.63731 -1.22396 0.405051 -0.54999 0.43997
## 8 0.5583 0.63731 -1.22396 -1.029299 1.29069 -0.66249
## 9 0.5583 0.63731 -1.22396 -0.400662 0.94009 -1.66619
## 10 0.4313 1.32739 1.70492 0.094693 -0.26410 0.77515
## 11 0.4313 1.32739 1.70492 0.132854 0.63174 -0.44504
## 12 0.5601 -1.38599 0.35086 0.207790 1.17765 -0.18459
## 13 0.5601 -1.38599 0.35086 0.002875 0.94375 -0.06831
## 14 -1.8785 -0.05503 -0.06852 -1.775732 -0.79413 1.57589
## 15 -1.8785 -0.05503 -0.06852 -1.716492 -0.23149 0.29543
## 16 0.5601 -1.38599 0.35086 -1.640220 0.82453 1.87405
## 17 -1.8785 -0.05503 -0.06852 2.662063 0.56436 2.53853
## 18 -1.8785 -0.05503 -0.06852 0.605071 -1.84599 -1.92760
## 19 -1.8785 -0.05503 -0.06852 2.309910 2.21260 -0.08448
## 20 -1.8785 -0.05503 -0.06852 -1.476716 -0.09132 -0.90420
##
##
## Biplot scores for constraining variables
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## ManagementHF 0.3751 0.42813 -0.82221 0 0 0
## ManagementNM -0.9989 -0.02926 -0.03644 0 0 0
## ManagementSF 0.3648 -0.90262 0.22849 0 0 0
##
##
## Centroids for factor constraints
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## ManagementBF 0.4313 1.32739 1.70492 0 0 0
## ManagementHF 0.5583 0.63731 -1.22396 0 0 0
## ManagementNM -1.8785 -0.05503 -0.06852 0 0 0
## ManagementSF 0.5601 -1.38599 0.35086 0 0 0# Permutation test on CCA
MCperm <- permutest(dune.ccaMN, permutations=999)
MCperm##
## Permutation test for cca under reduced model
##
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = dune ~ Management, data = dune.env.original)
## Permutation test for all constrained eigenvalues
## Df Inertia F Pr(>F)
## Model 3 0.60384 2.1307 0.002 **
## Residual 16 1.51143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Plot the ordination with only Management
ordiplot(dune.ccaMN, type = "points", main = "CCA only Management")
ordihull(dune.ccaMN,groups = dune.env.original$Management,draw = "polygon",col = "grey70", label = T)
##Q16
Question 16: Next, you want to investigate whether initial soil characteristics (A1 and moisture) also are important for differences in species composition between fields? To do this, do another CCA with A1 and Moisture as environmental variables and perform a significance test.
Explain the results of the permutations! Compare the results in exercises 16 and 17!
# CCA using A1 and moisture as the environmental variable
dune.ccaA1ms <- cca(dune, dune.env[,1:2])
summary(dune.ccaA1ms)##
## Call:
## cca(X = dune, Y = dune.env[, 1:2])
##
## Partitioning of scaled Chi-square:
## Inertia Proportion
## Total 2.1153 1.0000
## Constrained 0.5357 0.2533
## Unconstrained 1.5795 0.7467
##
## Eigenvalues, and their contribution to the scaled Chi-square
##
## Importance of components:
## CCA1 CCA2 CA1 CA2 CA3 CA4 CA5
## Eigenvalue 0.4265 0.10926 0.4044 0.2809 0.17168 0.14665 0.10004
## Proportion Explained 0.2016 0.05165 0.1912 0.1328 0.08116 0.06933 0.04729
## Cumulative Proportion 0.2016 0.25327 0.4444 0.5772 0.65841 0.72774 0.77503
## CA6 CA7 CA8 CA9 CA10 CA11 CA12
## Eigenvalue 0.09194 0.07595 0.07326 0.06055 0.04690 0.04554 0.02269
## Proportion Explained 0.04347 0.03590 0.03464 0.02862 0.02217 0.02153 0.01073
## Cumulative Proportion 0.81850 0.85440 0.88904 0.91766 0.93983 0.96136 0.97209
## CA13 CA14 CA15 CA16 CA17
## Eigenvalue 0.020748 0.015104 0.010950 0.007971 0.004264
## Proportion Explained 0.009809 0.007141 0.005177 0.003768 0.002016
## Cumulative Proportion 0.981898 0.989039 0.994216 0.997984 1.000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## CCA1 CCA2
## Eigenvalue 0.4265 0.1093
## Proportion Explained 0.7961 0.2039
## Cumulative Proportion 0.7961 1.0000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
##
##
## Species scores
##
## CCA1 CCA2 CA1 CA2 CA3 CA4
## Achimill -0.80522 -0.034144 0.38619 -0.039691 -0.22546 0.96607
## Agrostol 0.80865 -0.056511 -0.47652 0.122059 0.14051 -0.06689
## Airaprae 0.45857 0.929764 3.17582 -1.067783 0.42954 0.56399
## Alopgeni 0.48206 0.212415 -0.72292 -0.643671 -0.07254 -0.26939
## Anthodor -0.35975 0.070548 1.37653 -0.370765 -0.43073 0.67331
## Bellpere -0.68384 -0.258107 -0.18712 -0.007862 0.64294 0.10512
## Bromhord -0.73703 -0.006405 -0.18723 -0.075226 0.34288 0.73667
## Chenalbu 1.29869 0.096750 -0.82155 -1.823846 -0.64420 -0.20148
## Cirsarve -0.32362 0.081849 -0.78087 -0.495313 1.57388 -0.13821
## Comapalu 1.87488 -2.431666 0.31732 0.811112 0.65790 0.62535
## Eleopalu 1.40921 -0.388246 -0.16553 1.155390 -0.03457 0.28163
## Elymrepe -0.49671 0.018245 -0.62108 -0.338009 0.53782 0.33711
## Empenigr 0.99750 1.418422 3.18976 -1.141130 0.86374 -0.53939
## Hyporadi 0.23008 0.822734 2.58776 -0.729843 0.67116 -0.15454
## Juncarti 1.06605 0.324950 -0.35744 0.781955 -0.23317 -0.07351
## Juncbufo 0.38686 0.256698 -0.65114 -1.109108 -0.73373 -0.53684
## Lolipere -0.68296 0.121240 -0.23331 0.090110 0.25987 0.11064
## Planlanc -0.88359 -0.341888 0.54345 0.181775 -0.58711 -0.04526
## Poaprat -0.47409 0.147648 -0.23347 -0.057499 0.22428 0.06748
## Poatriv -0.12030 0.038928 -0.44519 -0.437623 -0.13930 0.16610
## Ranuflam 1.33049 -0.042806 -0.18729 0.945186 -0.15012 0.05862
## Rumeacet -0.64258 -0.500123 0.06444 -0.221997 -1.25444 -0.10685
## Sagiproc 0.34029 0.378747 -0.02196 -0.817623 0.34366 -0.60087
## Salirepe 0.45049 0.900374 1.06230 1.233516 0.23704 -1.16588
## Scorautu -0.09368 -0.002416 0.42515 -0.010154 0.17886 -0.20059
## Trifprat -0.98927 -0.564091 0.29087 0.151495 -1.58839 0.02149
## Trifrepe -0.03239 -0.209524 0.05308 -0.020468 0.01832 0.13533
## Vicilath -0.90579 -0.038667 0.46915 0.605435 0.77606 -1.05955
## Bracruta 0.05637 -0.129721 0.19033 0.331955 -0.10193 -0.58438
## Callcusp 1.36154 -0.179077 -0.05984 1.635569 -0.03460 0.51697
##
##
## Site scores (weighted averages of species scores)
##
## CCA1 CCA2 CA1 CA2 CA3 CA4
## 1 -1.2649 0.791185 -0.93676 0.24888 1.48190 1.10604
## 2 -0.8498 0.001097 -0.52316 -0.01984 1.10655 1.04020
## 3 -0.2389 0.419553 -0.88536 -0.49083 0.84021 -0.07468
## 4 -0.1614 0.505509 -0.78087 -0.49531 1.57388 -0.13821
## 5 -1.1319 -1.169583 0.38927 -0.40435 -0.83372 0.44207
## 6 -1.0755 -1.505312 0.37439 0.32710 -1.98034 -0.25573
## 7 -1.0659 -0.672190 -0.01631 0.26832 -1.36315 0.29399
## 8 0.7940 0.411933 -0.62913 0.12266 -0.24099 0.27183
## 9 0.1152 0.609288 -0.86673 -0.69040 -0.44865 -0.12278
## 10 -1.0563 -0.321545 0.22021 0.17724 0.09727 1.45794
## 11 -0.7576 0.451831 0.51559 0.52629 1.08266 -1.56601
## 12 0.5476 0.527784 -0.62514 -1.68048 -0.77125 -1.61783
## 13 0.7085 0.704490 -0.82155 -1.82385 -0.64420 -0.20148
## 14 1.9879 -3.320249 0.25773 1.07992 0.64829 1.44833
## 15 2.0283 -2.735712 0.37691 0.54230 0.66751 -0.19763
## 16 2.0458 -0.760204 -0.59118 1.32873 -0.47799 0.49495
## 17 -0.6410 1.939122 3.15491 -0.95776 -0.22177 2.21905
## 18 -0.5693 0.202225 0.62521 1.19191 0.84167 -2.56413
## 19 0.3779 3.849064 3.18976 -1.14113 0.86374 -0.53939
## 20 1.9395 0.807201 0.04808 2.68327 -0.50175 -0.70282
##
##
## Site constraints (linear combinations of constraining variables)
##
## CCA1 CCA2 CA1 CA2 CA3 CA4
## 1 -1.2002 0.36172 -0.93676 0.24888 1.48190 1.10604
## 2 -1.1086 -0.04053 -0.52316 -0.01984 1.10655 1.04020
## 3 -0.3105 0.02438 -0.88536 -0.49083 0.84021 -0.07468
## 4 -0.3236 0.08185 -0.78087 -0.49531 1.57388 -0.13821
## 5 -0.7419 -1.64952 0.38927 -0.40435 -0.83372 0.44207
## 6 -1.0038 -0.50024 0.37439 0.32710 -1.98034 -0.25573
## 7 -1.2002 0.36172 -0.01631 0.26832 -1.36315 0.29399
## 8 1.0630 1.13110 -0.62913 0.12266 -0.24099 0.27183
## 9 0.3042 0.89380 -0.86673 -0.69040 -0.44865 -0.12278
## 10 -0.4415 0.59902 0.22021 0.17724 0.09727 1.45794
## 11 -1.1086 -0.04053 0.51559 0.52629 1.08266 -1.56601
## 12 0.5792 -0.31295 -0.62514 -1.68048 -0.77125 -1.61783
## 13 1.2987 0.09675 -0.82155 -1.82385 -0.64420 -0.20148
## 14 1.7308 -1.79956 0.25773 1.07992 0.64829 1.44833
## 15 2.0189 -3.06377 0.37691 0.54230 0.66751 -0.19763
## 16 1.2594 0.26914 -0.59118 1.32873 -0.47799 0.49495
## 17 -0.3498 0.19678 3.15491 -0.95776 -0.22177 2.21905
## 18 -0.9645 -0.67263 0.62521 1.19191 0.84167 -2.56413
## 19 0.9975 1.41842 3.18976 -1.14113 0.86374 -0.53939
## 20 0.9713 1.53335 0.04808 2.68327 -0.50175 -0.70282
##
##
## Biplot scores for constraining variables
##
## CCA1 CCA2 CA1 CA2 CA3 CA4
## A1 0.6034 -0.7974 0 0 0 0
## Moisture 0.9750 0.2222 0 0 0 0# Permutation test on CCA
MCpermA1ms <- permutest(dune.ccaA1ms, permutations=999)
MCpermA1ms ##
## Permutation test for cca under reduced model
##
## Permutation: free
## Number of permutations: 999
##
## Model: cca(X = dune, Y = dune.env[, 1:2])
## Permutation test for all constrained eigenvalues
## Df Inertia F Pr(>F)
## Model 2 0.53573 2.8829 0.001 ***
## Residual 17 1.57954
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Plot the ordination
ordiplot(dune.ccaA1ms, type="text", main="CCA A1 and moisture")
Direct methods with testing of environmental variables
##Q17
Question 17: Partial ordination
Now, you have concluded that both soil characteristics and management type are determining the species composition significantly. But, what if the observed differences between management types are not caused by management type but by initial differences in soil characteristics? To investigate this, you have to test whether there still is a difference in vegetation between management types, after accounting for (i.e. removing) the effects of soil characteristics. This is done by a partial CCA ordination. Partial ordinations are used to eliminate effects of selected variables by specifying them as covariates. In this case, use the two soil characterising variables (A1 and Moisture) as covariates, and the four management types as environmental variables.
What does the permutation test tell you? Also, compare the results with the results from the permutation test on only management types!
#This partials out the effect of A1 and Moisture before analysing the effects of management
MN.cca <- cca(dune ~ Management + Condition(A1 + Moisture), data=dune.env.original)
summary(MN.cca)##
## Call:
## cca(formula = dune ~ Management + Condition(A1 + Moisture), data = dune.env.original)
##
## Partitioning of scaled Chi-square:
## Inertia Proportion
## Total 2.1153 1.0000
## Conditioned 0.7437 0.3516
## Constrained 0.3954 0.1869
## Unconstrained 0.9761 0.4615
##
## Eigenvalues, and their contribution to the scaled Chi-square
## after removing the contribution of conditiniong variables
##
## Importance of components:
## CCA1 CCA2 CCA3 CA1 CA2 CA3 CA4
## Eigenvalue 0.2385 0.08653 0.07036 0.3064 0.13191 0.11516 0.10947
## Proportion Explained 0.1739 0.06309 0.05130 0.2234 0.09618 0.08396 0.07982
## Cumulative Proportion 0.1739 0.23701 0.28831 0.5117 0.60787 0.69183 0.77165
## CA5 CA6 CA7 CA8 CA9 CA10 CA11
## Eigenvalue 0.07724 0.07575 0.04871 0.03758 0.03106 0.02102 0.012542
## Proportion Explained 0.05632 0.05523 0.03552 0.02740 0.02264 0.01533 0.009145
## Cumulative Proportion 0.82796 0.88320 0.91872 0.94612 0.96876 0.98409 0.993236
## CA12
## Eigenvalue 0.009277
## Proportion Explained 0.006764
## Cumulative Proportion 1.000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## CCA1 CCA2 CCA3
## Eigenvalue 0.2385 0.08653 0.07036
## Proportion Explained 0.6032 0.21882 0.17794
## Cumulative Proportion 0.6032 0.82206 1.00000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
##
##
## Species scores
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## Achimill -0.22184 0.1625650 -0.486532 -0.29895 -0.1910403 0.80559
## Agrostol 0.31801 -0.1591008 0.237934 0.39533 -0.0572108 -0.04906
## Airaprae -1.90497 0.0498869 0.657458 -2.65101 -0.0560655 0.69574
## Alopgeni 0.75017 -0.2904809 0.198070 -0.05494 -0.0179481 -0.14369
## Anthodor -0.67482 0.5664611 -0.023186 -1.13302 -0.4598130 0.53783
## Bellpere -0.11367 -0.2763144 -0.178950 0.22586 0.4032814 -0.05454
## Bromhord 0.11346 -0.0597708 -0.723087 0.06592 0.2219471 0.35933
## Chenalbu 1.65366 -0.8426420 0.255191 -1.21491 0.5419363 -0.08793
## Cirsarve 0.59767 -0.1296937 0.958252 0.37466 0.7316043 -0.84542
## Comapalu -0.99169 0.0882190 -0.518567 0.58438 0.4791343 0.38437
## Eleopalu 0.05742 -0.0080668 -0.116978 0.81272 -0.4722220 0.23229
## Elymrepe 0.29577 -0.0331848 0.246873 0.27977 1.0085693 0.45803
## Empenigr -1.55294 -0.2816661 0.504499 -2.90039 0.7122104 -0.36569
## Hyporadi -1.40452 -0.3342659 0.124456 -2.19156 -0.0002564 -0.14329
## Juncarti -0.11253 0.4647145 -0.179642 0.86164 0.0195413 -0.14119
## Juncbufo 0.43106 -0.1028278 0.177721 -0.30203 -0.0236828 0.02339
## Lolipere 0.21702 -0.0888942 -0.060884 0.16623 0.2453692 -0.06121
## Planlanc -0.28302 0.4634731 0.004689 -0.25092 -0.5717915 -0.09677
## Poaprat 0.13569 0.0008428 -0.009655 0.12030 0.3263968 -0.07370
## Poatriv 0.52883 0.0242477 -0.008592 -0.09088 0.2279744 0.14426
## Ranuflam -0.13726 -0.0519714 -0.010035 0.76566 -0.2142823 0.15268
## Rumeacet 0.26836 0.8823894 0.261872 -0.26470 -0.5038215 0.15287
## Sagiproc 0.22017 -0.2827029 0.171294 -0.54989 0.2153613 -0.62519
## Salirepe -1.69992 -0.4384927 0.770005 0.49821 -0.1961342 -0.48638
## Scorautu -0.39518 -0.0987393 -0.107125 -0.24614 0.0727362 -0.21337
## Trifprat 0.31106 1.1523633 0.368431 -0.26291 -0.8156335 0.12256
## Trifrepe -0.14521 -0.0264189 -0.408072 -0.03594 -0.0125782 0.09669
## Vicilath -0.60359 -0.8681108 -0.984615 0.09323 -0.4824554 -1.41297
## Bracruta -0.26165 -0.0801219 0.147628 0.14077 -0.3547347 -0.48658
## Callcusp -0.41589 -0.3345768 0.110571 1.18584 -0.8150813 1.08169
##
##
## Site scores (weighted averages of species scores)
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## 1 1.15369 -1.27332 0.4310 0.61137 2.18399 1.89365
## 2 0.61229 -1.61462 -1.6552 0.21908 1.29887 1.22585
## 3 0.91451 -0.84774 1.1194 0.39512 0.53538 -0.23449
## 4 0.74161 -0.94528 1.3644 0.37466 0.73160 -0.84542
## 5 -0.08086 1.40948 0.7508 -0.67144 0.18945 0.50159
## 6 0.02924 2.27467 0.8245 -0.24323 -1.32960 -0.14012
## 7 0.39168 1.23604 -0.2425 0.09643 -0.53580 0.40023
## 8 1.25965 0.57148 -0.9681 0.29444 0.79619 -1.07650
## 9 -0.05559 0.67863 -0.4734 0.59314 1.07760 0.29067
## 10 -0.48177 0.90282 -2.7160 0.03464 -0.84210 0.30477
## 11 -0.62527 -2.12839 -0.7924 -0.33410 -0.57320 -2.01847
## 12 0.06671 -0.81435 0.5681 -0.71177 -1.29312 -0.34880
## 13 1.98505 -0.59361 0.3088 -1.21491 0.54194 -0.08793
## 14 -0.61445 -0.66391 -1.6977 0.80540 0.13905 1.83177
## 15 -0.65898 0.29387 -0.4360 0.36335 0.81922 -1.06304
## 16 0.92510 0.06396 0.1286 0.64651 -2.00831 0.86205
## 17 -3.28246 2.50840 0.7073 -2.27695 -1.20848 2.28790
## 18 -1.48400 -2.23671 0.9243 1.00646 0.05868 -1.91972
## 19 -2.94792 0.05320 1.4111 -2.90039 0.71221 -0.36569
## 20 -0.81071 0.06920 1.0103 2.23242 -0.89403 0.30122
##
##
## Site constraints (linear combinations of constraining variables)
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## 1 0.799198 -1.62562 1.8039 0.61137 2.18399 1.89365
## 2 -0.004951 -1.34725 -1.5881 0.21908 1.29887 1.22585
## 3 0.606048 -0.12417 0.9430 0.39512 0.53538 -0.23449
## 4 0.597671 -0.12969 0.9583 0.37466 0.73160 -0.84542
## 5 0.469289 1.25664 0.0800 -0.67144 0.18945 0.50159
## 6 0.301753 1.14623 0.3854 -0.24323 -1.32960 -0.14012
## 7 0.176100 1.06342 0.6144 0.09643 -0.53580 0.40023
## 8 0.879780 1.74702 -0.6594 0.29444 0.79619 -1.07650
## 9 -0.363187 1.16959 -0.3949 0.59314 1.07760 0.29067
## 10 -0.340508 0.06035 -2.1894 0.03464 -0.84210 0.30477
## 11 -0.004951 -1.34725 -1.5881 -0.33410 -0.57320 -2.01847
## 12 0.435824 -1.40351 0.4739 -0.71177 -1.29312 -0.34880
## 13 1.653661 -0.84264 0.2552 -1.21491 0.54194 -0.08793
## 14 -1.083840 0.02749 -0.3506 0.80540 0.13905 1.83177
## 15 -0.899549 0.14895 -0.6865 0.36335 0.81922 -1.06304
## 16 1.628531 -0.85920 0.3010 0.64651 -2.00831 0.86205
## 17 -2.433019 0.54722 0.8869 -2.27695 -1.20848 2.28790
## 18 -2.063955 -0.83829 1.4271 1.00646 0.05868 -1.91972
## 19 -1.552942 -0.28167 0.5045 -2.90039 0.71221 -0.36569
## 20 -1.569696 -0.29271 0.5350 2.23242 -0.89403 0.30122
##
##
## Biplot scores for constraining variables
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## ManagementHF 0.1944 0.85341 0.01122 0 0 0
## ManagementNM -0.8301 -0.09912 0.21079 0 0 0
## ManagementSF 0.6036 -0.46400 0.47327 0 0 0
##
##
## Centroids for factor constraints
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## ManagementBF -0.1283 -0.8299 -1.8091 0 0 0
## ManagementHF 0.2894 1.2704 0.0167 0 0 0
## ManagementNM -1.5611 -0.1864 0.3964 0 0 0
## ManagementSF 0.9269 -0.7125 0.7267 0 0 0#permutation test
permutest(MN.cca, permutations=999)##
## Permutation test for cca under reduced model
##
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = dune ~ Management + Condition(A1 + Moisture), data
## = dune.env.original)
## Permutation test for all constrained eigenvalues
## Df Inertia F Pr(>F)
## Model 3 0.39543 1.6205 0.021 *
## Residual 12 0.97610
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Plot the ordination
ordiplot(MN.cca, type = "points")
ordihull(MN.cca,groups = dune.env.original$Management,draw = "polygon",col = "grey70", label = T)
##Q18
Question 18: The amount of applied manure is a third important variable given in the environmental data set associated with the dune meadow data set. Repeat exercise 17, but with manure added to the set of covariates.
Inspect and interpret the results of the permutations!
#This partials out the effect of A1, Moisture, and Manure before analysing the effects of management
MN.cca2 <- cca(dune ~ Management + Condition(A1 + Moisture + Manure), data=dune.env.original)
summary(MN.cca2)##
## Call:
## cca(formula = dune ~ Management + Condition(A1 + Moisture + Manure), data = dune.env.original)
##
## Partitioning of scaled Chi-square:
## Inertia Proportion
## Total 2.1153 1.00000
## Conditioned 1.2207 0.57708
## Constrained 0.1524 0.07204
## Unconstrained 0.7422 0.35088
##
## Eigenvalues, and their contribution to the scaled Chi-square
## after removing the contribution of conditiniong variables
##
## Importance of components:
## CCA1 CCA2 CA1 CA2 CA3 CA4 CA5
## Eigenvalue 0.1019 0.05047 0.2759 0.1130 0.1047 0.07949 0.06436
## Proportion Explained 0.1139 0.05642 0.3084 0.1263 0.1170 0.08885 0.07195
## Cumulative Proportion 0.1139 0.17035 0.4787 0.6051 0.7221 0.81095 0.88290
## CA6 CA7 CA8 CA9
## Eigenvalue 0.04147 0.03105 0.01842 0.01381
## Proportion Explained 0.04636 0.03471 0.02059 0.01543
## Cumulative Proportion 0.92926 0.96397 0.98457 1.00000
##
## Accumulated constrained eigenvalues
## Importance of components:
## CCA1 CCA2
## Eigenvalue 0.1019 0.05047
## Proportion Explained 0.6688 0.33119
## Cumulative Proportion 0.6688 1.00000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
##
##
## Species scores
##
## CCA1 CCA2 CA1 CA2 CA3 CA4
## Achimill -0.133335 0.120821 -0.11336 0.6510292 -0.337080 0.22393
## Agrostol -0.059233 -0.132637 0.37444 -0.0370997 0.095375 -0.07786
## Airaprae 0.407915 -0.230470 -2.44958 1.0332387 0.418358 0.36362
## Alopgeni -0.280823 -0.253038 -0.01362 -0.2317652 0.041553 0.10891
## Anthodor 0.639819 -0.006334 -0.83017 0.5633922 -0.000667 0.16112
## Bellpere -0.434577 0.069176 0.21427 -0.3756168 -0.309418 0.31025
## Bromhord -0.599725 0.321928 0.10185 0.0261899 -0.332718 -0.04073
## Chenalbu 0.057009 -1.809598 -0.99123 -1.3005988 -2.653184 0.02770
## Cirsarve 0.055897 0.173517 0.05198 -0.8024478 0.546950 -0.82929
## Comapalu -0.596782 0.942411 0.16931 0.5922871 -0.074486 -0.77810
## Eleopalu -0.128649 0.070611 0.82317 0.5764876 0.621879 0.16780
## Elymrepe -0.006479 0.312140 -0.04221 0.3913889 -0.497683 0.46659
## Empenigr 0.356332 -0.139879 -2.96848 -0.1566239 0.522892 -0.37091
## Hyporadi 0.131937 -0.381471 -2.18339 0.4041744 0.743463 -0.02112
## Juncarti 0.168285 0.541592 0.75315 0.0785039 0.371408 0.73033
## Juncbufo -0.059236 -0.346550 -0.24256 -0.2636189 -0.567659 -0.18092
## Lolipere -0.124276 0.092667 -0.01431 0.0985747 0.046798 -0.02100
## Planlanc 0.515751 -0.077257 -0.06547 -0.0155724 0.129069 -0.06399
## Poaprat -0.080901 0.116929 -0.04072 -0.0256861 -0.125095 0.11124
## Poatriv -0.056841 -0.011211 -0.08430 -0.1019155 -0.424954 0.07988
## Ranuflam -0.039091 0.005367 0.81474 0.0492154 -0.156240 0.03827
## Rumeacet 0.974169 0.193476 -0.00166 0.0550963 0.027334 -0.19856
## Sagiproc -0.127410 -0.181867 -0.63145 -0.5495718 0.316077 -0.35896
## Salirepe 0.248228 -0.408801 0.66350 -0.7497857 -0.070381 0.23967
## Scorautu -0.144092 0.019065 -0.26341 -0.2057033 0.063209 -0.02494
## Trifprat 1.300548 0.240744 0.09341 0.0640604 0.167503 -0.48798
## Trifrepe -0.190514 0.114486 -0.02230 0.0004746 -0.178668 -0.46336
## Vicilath -0.479706 -0.825859 -0.01391 -0.8227623 0.713571 -0.16779
## Bracruta 0.143858 -0.215431 0.20953 -0.3155810 0.462924 0.10961
## Callcusp -0.030780 -0.303185 1.35836 1.1001438 -0.074428 -1.09279
##
##
## Site scores (weighted averages of species scores)
##
## CCA1 CCA2 CA1 CA2 CA3 CA4
## 1 0.076849 1.6054 -0.38254 2.52074 -1.15972 1.05976
## 2 -2.949104 1.3252 0.12543 0.38912 -0.76331 0.55160
## 3 0.000495 -0.5491 0.11367 -0.23158 -0.09345 0.45606
## 4 -0.031180 -0.1541 0.05198 -0.80245 0.54695 -0.82929
## 5 1.275980 0.1262 -0.37104 0.08451 -0.62045 0.96726
## 6 1.992850 0.1218 0.33239 -0.07570 0.55582 -0.86651
## 7 0.458162 0.9887 -0.03962 0.39302 -0.01533 -0.99690
## 8 -0.341512 1.4118 -0.09208 -0.80160 0.81681 0.41771
## 9 0.634800 1.5936 0.12543 0.38912 -0.76331 0.55160
## 10 -0.403563 0.6018 0.42285 0.07216 -0.57670 -0.06757
## 11 -0.290887 -2.9002 -0.73283 -0.60769 1.77679 -0.63318
## 12 -0.761760 -1.9123 -0.15051 -0.46695 0.91598 -0.66193
## 13 -0.306982 -2.4614 -0.99123 -1.30060 -2.65318 0.02770
## 14 -0.984376 1.0251 0.54638 1.60362 -1.15108 -3.50187
## 15 0.113636 1.6486 -0.20777 -0.41904 1.00211 1.94566
## 16 0.165589 -0.4483 1.15086 1.79585 1.68170 0.67435
## 17 1.249098 0.2014 -1.67122 2.81803 0.26156 1.46543
## 18 -1.372702 -1.5766 0.98715 -2.14783 -0.12260 0.66276
## 19 1.211626 -1.3634 -2.96848 -0.15662 0.52289 -0.37091
## 20 0.057337 0.6223 2.64850 -0.26686 -0.39501 0.35216
##
##
## Site constraints (linear combinations of constraining variables)
##
## CCA1 CCA2 CA1 CA2 CA3 CA4
## 1 -0.232348 -0.85528 -0.38254 2.52074 -1.15972 1.05976
## 2 -2.714892 0.72966 0.12543 0.38912 -0.76331 0.55160
## 3 0.041672 0.18967 0.11367 -0.23158 -0.09345 0.45606
## 4 0.055897 0.17352 0.05198 -0.80245 0.54695 -0.82929
## 5 1.504037 0.29428 -0.37104 0.08451 -0.62045 0.96726
## 6 1.788549 -0.02879 0.33239 -0.07570 0.55582 -0.86651
## 7 -0.122942 0.86104 -0.03962 0.39302 -0.01533 -0.99690
## 8 -0.006331 2.16478 -0.09208 -0.80160 0.81681 0.41771
## 9 0.869012 0.99804 0.12543 0.38912 -0.76331 0.55160
## 10 -0.266549 -0.23023 0.42285 0.07216 -0.57670 -0.06757
## 11 -0.782403 -1.00056 -0.73283 -0.60769 1.77679 -0.63318
## 12 -1.042814 -1.19765 -0.15051 -0.46695 0.91598 -0.66193
## 13 0.057009 -1.80960 -0.99123 -1.30060 -2.65318 0.02770
## 14 -0.440301 0.76472 0.54638 1.60362 -1.15108 -3.50187
## 15 -0.753264 1.12010 -0.20777 -0.41904 1.00211 1.94566
## 16 0.099686 -1.85806 1.15086 1.79585 1.68170 0.67435
## 17 0.485289 -0.36636 -1.67122 2.81803 0.26156 1.46543
## 18 -0.087468 -1.07208 0.98715 -2.14783 -0.12260 0.66276
## 19 0.356332 -0.13988 -2.96848 -0.15662 0.52289 -0.37091
## 20 0.384783 -0.17219 2.64850 -0.26686 -0.39501 0.35216
##
##
## Biplot scores for constraining variables
##
## CCA1 CCA2 CA1 CA2 CA3 CA4
## ManagementHF 0.5735 0.5495 0 0 0 0
## ManagementSF -0.1000 -0.5202 0 0 0 0
##
##
## Centroids for factor constraints
##
## CCA1 CCA2 CA1 CA2 CA3 CA4
## ManagementBF -1.2865 -0.09634 0 0 0 0
## ManagementHF 0.8538 0.81795 0 0 0 0
## ManagementSF -0.1536 -0.79878 0 0 0 0#permutation test
permutest(MN.cca2, permutations=999)##
## Permutation test for cca under reduced model
##
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = dune ~ Management + Condition(A1 + Moisture +
## Manure), data = dune.env.original)
## Permutation test for all constrained eigenvalues
## Df Inertia F Pr(>F)
## Model 2 0.15239 0.9239 0.542
## Residual 9 0.74220
##Q19
Question 19: Forward selection
We have now concluded that management type is of subordinate importance. Now, you want to test the importance or significance of individual environmental variables. This is done by Forward selection. What is the effect of Forward selection? In what way do these results differ from the results obtained when the whole Environmental data table was used as explanatory variables? Were the set of explanatory variables selected by the Forward selection the best set? Why would one perform a p-value correction?
# A first CCA with all environmental variables
cca.dune.all <- cca(dune ~ ., data = dune.env) # the dot after the tilde (~) means "include all from data"
anova(cca.dune.all) # Test overall significance of the explanatory variables. If not significant no need to proceed with step wise selectionadjR2.cca.dune <- RsquareAdj(cca.dune.all)$adj.r.squared # How much the adjusted R2 explains when all variables in dune.env are included
adjR2.cca.dune## [1] 0.2299496# A null model with only intercept, needed for the step wise selection
cca.dune.0 <- cca(dune ~ 1, data = dune.env) # model containing only species matrix and intercept
# The step wise selection, based on R2
sel.env.dune <- ordiR2step(cca.dune.0, scope = formula (cca.dune.all), R2scope = adjR2.cca.dune, direction = 'forward', permutations = 999)## Step: R2.adj= 0
## Call: dune ~ 1
##
## R2.adjusted
## <All variables> 0.23030613
## + Management 0.15184053
## + Moisture 0.15084705
## + Manure 0.10001447
## + A1 0.05470899
## + Use 0.01489235
## <none> 0.00000000
##
## Df AIC F Pr(>F)
## + Management 3 86.935 2.1307 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.1523214
## Call: dune ~ Management
##
## R2.adjusted
## + Moisture 0.2462733
## <All variables> 0.2303061
## + A1 0.2000052
## + Use 0.1532069
## <none> 0.1523214
## + Manure 0.1466277sel.env.dune## Call: cca(formula = dune ~ Management, data = dune.env)
##
## Inertia Proportion Rank
## Total 2.1153 1.0000
## Constrained 0.6038 0.2855 3
## Unconstrained 1.5114 0.7145 16
## Inertia is scaled Chi-square
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3
## 0.3186 0.1825 0.1027
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 CA9 CA10 CA11
## 0.4474 0.2030 0.1630 0.1346 0.1294 0.0949 0.0790 0.0653 0.0500 0.0432 0.0387
## CA12 CA13 CA14 CA15 CA16
## 0.0239 0.0177 0.0092 0.0080 0.0042sel.env.dune$anovaSince there is (a potentially high) number of tests of significance during the forward selection procedure, it is better to apply a correction for multiple testing issue. Here we apply the Holm correction.
sel.env.dune_adj <- sel.env.dune
sel.env.dune_adj$anova$`Pr(>F)` <- p.adjust(sel.env.dune$anova$`Pr(>F)`, method = 'holm', n = ncol(dune.env)) # 7 other methods available
sel.env.dune_adj$anova
##Q20
Question 20: Interpreting ordination results using species traits
Experienced plant ecologists may already have looked at the species in the Dune Meadow ordination graphs and concluded that species with similar traits occur together. This is possible if you have good knowledge about plant species ecological preferences, and if there are relatively few species in your dataset. In this exercise we will use tabulated data on species ecological preferences (the Ellenberg indicator values) to interpret the results of the ordinations. Which Ellenberg values are most important for the distribution of the species? What do the different axes represent in terms of environmental gradients?
#load Ellenberg data
dune.ell <- readRDS("Ellenberg.RDS")
dune.mean.ell <- readRDS("Mean_Ellenberg.RDS")
#The vegan implementation of forward selection of variables does not accept NA values.
#To get around this we will replace the two NAs with the mean value for all plots in that column.
#Do not do this in a real analysis without seriously thinking about how it affects your results!
## Code to replace missing values with the column mean, and save the results as a new dataframe
## called dune.mean.ell.impute (since we are imputing the missing values).
dune.mean.ell.impute <- data.frame(
sapply(
dune.mean.ell,
function(x) ifelse(is.na(x),
mean(x, na.rm = TRUE),
x)))
#CCA analysis
#create global model with CCA (including all variables) and test it's significance, and if it is significant,
#we use the ordistep function with appropriate arguments to do forward selection of variables.
cca1 <- cca(dune ~ ., data = dune.mean.ell.impute) # full model (with all explanatory variables)
anova(cca1) #overall model is significant# Start by making ordinations
cca0 <- cca (dune~ 1, data = dune.mean.ell.impute) # empty model only with intercept
cca1 <- cca(dune ~ ., data = dune.mean.ell.impute, na.action = na.omit) # full model (with all explanatory variables)
cca1## Call: cca(formula = dune ~ L + T + K + F + R + N + S, data =
## dune.mean.ell.impute, na.action = na.omit)
##
## Inertia Proportion Rank
## Total 2.1153 1.0000
## Constrained 1.3313 0.6294 7
## Unconstrained 0.7840 0.3706 12
## Inertia is scaled Chi-square
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7
## 0.4646 0.3298 0.1774 0.1490 0.1091 0.0597 0.0417
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 CA9 CA10
## 0.24495 0.14477 0.08921 0.07974 0.06167 0.04827 0.03879 0.02701 0.01952 0.01271
## CA11 CA12
## 0.01006 0.00729plot(cca1, main="CCA, all data")
#Stepwise approach, using "ordistep"
step2<-ordistep(cca0, scope = formula(cca1), direction="forward")##
## Start: dune ~ 1
##
## Df AIC F Pr(>F)
## + N 1 85.551 4.1029 0.005 **
## + L 1 86.138 3.4630 0.005 **
## + F 1 86.635 2.9367 0.005 **
## + T 1 86.624 2.9479 0.010 **
## + S 1 87.111 2.4440 0.015 *
## + R 1 87.551 1.9994 0.025 *
## + K 1 87.850 1.7028 0.080 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: dune ~ N
##
## Df AIC F Pr(>F)
## + F 1 83.748 3.5595 0.005 **
## + S 1 84.490 2.8114 0.005 **
## + L 1 84.501 2.8000 0.005 **
## + R 1 84.753 2.5526 0.015 *
## + T 1 85.234 2.0874 0.020 *
## + K 1 85.992 1.3778 0.160
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: dune ~ N + F
##
## Df AIC F Pr(>F)
## + T 1 82.828 2.5154 0.020 *
## + R 1 83.824 1.6159 0.060 .
## + S 1 83.715 1.7119 0.075 .
## + L 1 83.543 1.8649 0.080 .
## + K 1 83.881 1.5660 0.100 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: dune ~ N + F + T
##
## Df AIC F Pr(>F)
## + L 1 82.381 1.9525 0.015 *
## + S 1 82.728 1.6611 0.065 .
## + R 1 82.582 1.7830 0.070 .
## + K 1 83.052 1.3932 0.170
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: dune ~ N + F + T + L
##
## Df AIC F Pr(>F)
## + S 1 81.835 1.9003 0.055 .
## + R 1 82.299 1.5355 0.120
## + K 1 82.328 1.5135 0.145
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1step2$anova # which three are most significant in the final selection?##Q21
Question 21: Decomposition of variance
Perform the analyses and interpret the results! How much of the total variation (% of All) is explained by: - uniquely by Management (effect of Soil removed), - uniquely by Soil (effect of Management removed), - jointly by Soil and Management (the interaction) - Not by Soil or Management?
In usual analysis of variance experiments, the variance is decomposed into components. The same can be done in multivariate analyses. In this exercise you will decompose the variance in the Dune meadow data set into different variance components. You will use two groups of variance components: Group 1 is Management and Group 2 is the soil variables (A1 and Moisture). In ordinations, the variance is expressed by the sum of the eigenvalues.
The result gives you the fraction of the total variation that is explained by: a) uniquely by Management (effect of soil removed), b) uniquely by soil (effect of Management removed), c) jointly by Soil and Management (the interaction).
The total explained variation is the sum of a), b) and c).
The total variation in the dataset is the sum of all unconstrained eigenvalues.
data("dune")
data("dune.env")
#as we did in Question 3:PCA, we need to reformulate the environmental data to make
#new columns for each of the management methods.
dummy_management <- as.data.frame(model.matrix( ~ Management - 1, data=dune.env ))
#add these to the dataset
dune.env <- dune.env %>% select(A1, Moisture, Manure, Use) %>% cbind(.,dummy_management)
dune.env$Moisture <- as.numeric(as.character(dune.env$Moisture)) #make numeric
dune.env$Manure <- as.numeric(as.character(dune.env$Manure))
dune.env$Use <- as.numeric(dune.env$Use)
#make column names shorter
dune.env <- dune.env %>% rename(BF = ManagementBF, HF = ManagementHF,
NM = ManagementNM, SF = ManagementSF)
## variance partitioning
management <- dune.env[,c("BF", "HF", "NM", "SF")]
soil <- dune.env[,c("A1", "Moisture")]
# examine the explanatory variable of each class of variables.
varp <- varpart(dune, management, soil)
varp##
## Partition of variance in RDA
##
## Call: varpart(Y = dune, X = management, soil)
##
## Explanatory tables:
## X1: management
## X2: soil
##
## No. of explanatory tables: 2
## Total variation (SS): 1598.4
## Variance: 84.124
## No. of observations: 20
##
## Partition table:
## Df R.squared Adj.R.squared Testable
## [a+b] = X1 3 0.34747 0.22512 TRUE
## [b+c] = X2 2 0.25771 0.17038 TRUE
## [a+b+c] = X1+X2 5 0.51552 0.34249 TRUE
## Individual fractions
## [a] = X1|X2 3 0.17211 TRUE
## [b] 0 0.05301 FALSE
## [c] = X2|X1 2 0.11737 TRUE
## [d] = Residuals 0.65751 FALSE
## ---
## Use function 'rda' to test significance of fractions of interestplot (varp, digits = 2, Xnames = c('Management', 'Soil'), bg = c('red', 'blue'))
##Q22
Question 22: Analysing a time series with vegetation data
Is there a significant change in species composition over time? What is the effect of defining Plots as cofactors?
Quite often vegetation ecologists have data from repeated inventories in permanent plots. The overall question to answer is if there has been a significant and directional change in the species composition. In this exercise we have rearranged the Dune Meadow data so that it consists of only 10 plots, each analysed 2 times (same species as in all other exercises, just grouping and rearrangement of plots). We want to analyse if there has been a consistent change in species composition between the two inventories. We are thus not interested in differences between the ten plots.
#This partials out the effect of Plot before analysing the effects of Time
#First load the time series versions of the data
SpeTS <- readRDS("SpeTS.RDS")
EnvTS <- readRDS("EnvTS.RDS")
time.cca <- cca(SpeTS ~ Time + Condition(Plot), data=EnvTS)
time.cca## Call: cca(formula = SpeTS ~ Time + Condition(Plot), data = EnvTS)
##
## Inertia Proportion Rank
## Total 2.11526 1.00000
## Conditional 0.20908 0.09885 1
## Constrained 0.15941 0.07536 1
## Unconstrained 1.74677 0.82579 17
## Inertia is scaled Chi-square
##
## Eigenvalues for constrained axes:
## CCA1
## 0.15941
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.4215 0.2761 0.2461 0.1722 0.1358 0.0999 0.0913 0.0783
## (Showing 8 of 17 unconstrained eigenvalues)treat <- EnvTS$Time #
colvec <- c("red2", "green4") #set colours to be applied to different levels of factor "Time"
plot(time.cca, type = "n", display = c("lc"))#plots the axes
with(time.cca, points(time.cca, display = c("lc"), col = colvec[treat],#plots the points
pch = 21, xlim = c(-2,2), bg = colvec[treat]))
#permutation test
with(EnvTS, anova(time.cca, by="term", perm=500, strata=Plot))## Set of permutations < 'minperm'. Generating entire set.#Compare with results for time as only explanatory variable and no cofactors
time2.cca <- cca(SpeTS ~ Time, data=EnvTS)
time2.cca## Call: cca(formula = SpeTS ~ Time, data = EnvTS)
##
## Inertia Proportion Rank
## Total 2.11526 1.00000
## Constrained 0.14240 0.06732 1
## Unconstrained 1.97287 0.93268 18
## Inertia is scaled Chi-square
##
## Eigenvalues for constrained axes:
## CCA1
## 0.1424
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.4796 0.3721 0.2530 0.1736 0.1365 0.1063 0.0921 0.0783
## (Showing 8 of 18 unconstrained eigenvalues)#permutation test
with(EnvTS, anova(time2.cca, by="term", perm=500))## plot ellipsoid hulls
treat <- EnvTS$Time
plot(time2.cca, type = "n")#plots the axes
with(time2.cca, points(time2.cca, col = colvec[treat],#plots the points
pch = 21, bg = colvec[treat]))
ordihull(time2.cca,groups=treat,draw="polygon",col="grey70",label=T)#draw hulls around area of each level of factor "Time"
##Q23
Question 23: A multivariate Before-After-Control- Impact (BACI) study
Did the treatment have an effect? What was permuted?
In this exercise we will continue to use the rearranged Dune Meadow data from exercise 24. A difference is that the plots are now divided into four groups: 1. Control plots before a treatment (i.e. not treated) 2. Control plots after a treatment (i.e. not treated) 3. Impact plots before treatment (i.e. not treated) 4. Impact plots after treatment (this is where the treatment is made)
The four groups are indicated in variable Treat in the environmental data set (EnvTS). The question you want to answer is if the treatment caused a change in species composition that is significantly different from the change in the control plots.
baci.cca <- cca(SpeTS ~ Treat + Condition(Plot + Time), data=EnvTS)
baci.cca## Call: cca(formula = SpeTS ~ Treat + Condition(Plot + Time), data =
## EnvTS)
##
## Inertia Proportion Rank
## Total 2.1153 1.0000
## Conditional 0.3685 0.1742 2
## Constrained 0.2443 0.1155 1
## Unconstrained 1.5024 0.7103 16
## Inertia is scaled Chi-square
##
## Eigenvalues for constrained axes:
## CCA1
## 0.24432
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 CA9 CA10 CA11
## 0.3925 0.2472 0.1740 0.1453 0.1182 0.0933 0.0829 0.0570 0.0532 0.0441 0.0324
## CA12 CA13 CA14 CA15 CA16
## 0.0218 0.0141 0.0115 0.0082 0.0069with(EnvTS, anova(baci.cca, by="term", perm=500, strata=Treat))#Here with() is a special function that makes variables in dune.env visible to
#the following command. If you only type Moisture in an R prompt, you will get
#an error of missing variables## plot results
treat <- EnvTS$Treat
ordiplot(baci.cca,type="points")
ordihull(baci.cca,groups=treat,draw="polygon",col="grey70",label=T)#draw hulls around area of each level of factor "Time"
##Q24
Question 24: Procustes rotation and Co-correspondence analysis
Are the two data sets (species and bugs) correlated? Give a very short ecological interpretation of the results!
These are two ways of comparing ordinations or the structure in related multivariate datasets. They are used to compare how objects differ when they are described by different sets of descriptors, or if there are repeated samplings of the same objects using the same descriptors.
library(cocorresp)
# DCA on the plant species data
DCASpecies <- decorana(dune)
# DCA on the insect data (called Bugs)
Bugs <- readRDS("Bugs.RDS")
DCABugs <- decorana(Bugs)
DCABugs##
## Call:
## decorana(veg = Bugs)
##
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
##
## DCA1 DCA2 DCA3 DCA4
## Eigenvalues 0.6785 0.3236 0.17269 0.1292
## Decorana values 0.6853 0.2470 0.03452 0.0207
## Axis lengths 4.2552 2.2877 1.63208 1.6331# procrustes analysis on the two DCAs from different datasets
ProC <- procrustes(DCASpecies, DCABugs)
plot(ProC, kind=1, type="text")
# Labels show the position of the samples in the second ordination, and arrows point to their positions in the target ordination
# Function protest tests the non-randomness (`significance') between two configurations
ProCsig <- protest(dune.dca, DCABugs)
ProCsig##
## Call:
## protest(X = dune.dca, Y = DCABugs)
##
## Procrustes Sum of Squares (m12 squared): 0.4492
## Correlation in a symmetric Procrustes rotation: 0.7422
## Significance: 0.001
##
## Permutation: free
## Number of permutations: 999# Co-correspondence analysis
CoCor <- coca(dune ~ ., data = Bugs, method = "symmetric")
summary(CoCor)##
## Symmetric Co-Correspondence Analysis
##
## Call: symcoca(y = y, x = x, n.axes = n.axes, R0 = weights, symmetric =
## symmetric, nam.dat = nam.dat)
##
## Inertia:
## Total Explained Residual
## dune: 2.115 1.672 0.444
## Bugs: 2.358 2.280 0.078
##
## Eigenvalues:
## COCA 1 COCA 2 COCA 3 COCA 4 COCA 5 COCA 6 COCA 7
## 0.2751016 0.2551218 0.0662534 0.0210179 0.0075231 0.0035365 0.0014639
## COCA 8 COCA 9
## 0.0008311 0.0002375corAxis(CoCor)## COCA 1 COCA 2 COCA 3 COCA 4 COCA 5 COCA 6 COCA 7 COCA 8
## 0.9021743 0.9103017 0.8259081 0.8133787 0.8165470 0.7043982 0.7022924 0.5065464
## COCA 9
## 0.7012526##Q25
Question 25 : Classification
How do the different clustering techniques differ?
Start by testing the non-hierarchical K-means clustering, using Multivar / K-Means. Note that it is recommended that K-means need at least 100 observations (500 according to some sources) to be reliable! Let us ignore the sample size issue for the moment, and ask for 4 clusters (for later comparison with the four management types). Test different hierarchical agglomeration algorithms and similarity indices. Use at least the Euclidian and the Bray-Curtis similarity measures for the hierarchical clustering technique. Try also the Ward’s method!
library(dendextend)
#make a distance matrix using the Bray-Curtis method
dune.dist <- vegdist(dune, method = "bray")
# K-Means Cluster Analysis
fit <- kmeans(dune, 4, nstart = 50) # we specify a 4 cluster solution
#As the final result of k-means clustering result is sensitive to the random starting assignments, we specify nstart = 50. This means #that R will try 50 different random starting assignments and then select the best results corresponding to the one with the lowest #within cluster variation.
#Note that k-means clustering involves randomness so you won't neccessarily get exactly the same results if you repeat it.
#There are also several different algorithms that can be used in k-means clustering, and different software
#have can use different defaults.Check out the help file for the kmeans function if you want to know more about these.
#cluster number for each data point
fit$cluster## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 3 3 2 2 3 3 3 2 2 3 1 2 2 4 4 4 1 1 1 4#cluster sizes
fit$size## [1] 4 6 6 4#compare clusters to management category
table(fit$cluster, dune.env.original$Management)##
## BF HF NM SF
## 1 1 0 3 0
## 2 0 2 0 4
## 3 2 3 0 1
## 4 0 0 3 1#try some approaches to hierarchical clustering
dend2 <- dune %>% # data
dist(method = "euclidean") %>% # calculate a distance matrix, choose method
hclust(method = "ward.D") %>% # Hierarchical clustering, choose method
as.dendrogram # Turn the object into a dendrogram.
dend2 %>% set("labels_color") %>% set("branches_lwd") %>% set("branches_k_color", k = 4) %>% plot
#the dist() function does not include Bray-Curtis but we have already made a distance matrix using
#the vegdist() function from vegan.
dend3 <- dune.dist %>% # data (Bray-Curtis)
hclust(method = "aver") %>% # Hierarchical clustering, choose method
as.dendrogram # Turn the object into a dendrogram.
dend3 %>% set("labels_color") %>% set("branches_lwd") %>% set("branches_k_color", k = 4) %>% plot
#we can colour the labels based on their grouping in the kmeans analysis to compare clusters
dend3 %>% set("labels_color", fit$cluster) %>% set("branches_lwd") %>% set("branches_k_color", k = 4) %>% plot
#tanglegram is a nice function to compare dendrograms that you might want to try out as an extra
tanglegram(dend2, dend3)
#Q26 ### Question 26: Classification 2
Classification 2
#to follow##Q27
Question 27: Principal Response Curves
PRC
#to follow##Q28
Question 28: ANOSIM & SIMPER
What does the result tell you? What other types of predefined clusters could you use?
ANOSIM (ANalysis Of Similarities) is a non-parametric test of significant difference between two or more groups, based on any distance measure. In this case, use the clusters from the K-means exercise. Change to Bray-Curtis similarity index.
dune.dist <- vegdist(dune, method = "bray") #create distance matrix based on Bray-Curtis method
dune.ano <- anosim(dune.dist, fit$cluster) #compare groupings based on cluster analysis
summary(dune.ano)##
## Call:
## anosim(x = dune.dist, grouping = fit$cluster)
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.8226
## Significance: 0.001
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.116 0.171 0.208 0.263
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 12.5 76.875 115.50 153.50 188.0 148
## 1 10.0 62.250 66.25 101.75 134.0 6
## 2 3.0 14.250 22.00 29.50 53.5 15
## 3 1.0 6.000 23.50 55.75 91.0 15
## 4 6.0 15.250 19.50 32.00 60.0 6#plot(dune.ano)
#title(sub = "Anosim, groups = cluster analysis")
#ANOSIM gives you the P value and a R value. R value close to 1 indicates high separation between levels of your factor while R value #close to 0 indicate no separation between levels of your factor.dune.ano2 <- anosim(dune.dist, dune.env.original$Management) #try comparing groupings based on management instead
summary(dune.ano2)##
## Call:
## anosim(x = dune.dist, grouping = dune.env.original$Management)
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.2579
## Significance: 0.008
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.122 0.153 0.198 0.230
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 4 58.50 104.00 145.500 188.0 147
## BF 5 15.25 25.50 41.250 57.0 3
## HF 1 7.25 46.25 68.125 89.5 10
## NM 6 64.75 124.50 156.250 181.0 15
## SF 3 32.75 53.50 99.250 184.0 15#plot(dune.ano2, sub = "Anosim, groups = management")SIMPER
Compare the different groups and check which species that contributes mostly to the difference between pairs of clusters.
After a significant ANOSIM, you may want to know which species are primarily responsible for the observed difference between clusters. SIMPER (Similarity Percentage) will do this for you. The test does not come with significance testing. In the output tables, taxa are sorted in descending order of contribution to group difference.
sim <- simper(dune, fit$cluster) #start with kmeans clusters
summary(sim)##
## Contrast: 3_2
##
## average sd ratio ava avb cumsum
## Alopgeni 0.062264 0.033663 1.8496 0.3333 5.0000 0.1062
## Agrostol 0.060675 0.021567 2.8133 0.0000 4.6667 0.2097
## Lolipere 0.042051 0.035413 1.1875 5.3333 2.8333 0.2814
## Planlanc 0.036403 0.027638 1.3171 3.0000 0.0000 0.3435
## Sagiproc 0.032522 0.021454 1.5159 0.0000 2.5000 0.3990
## Elymrepe 0.030536 0.028642 1.0662 2.0000 2.3333 0.4511
## Achimill 0.029361 0.010739 2.7340 2.3333 0.0000 0.5012
## Rumeacet 0.028649 0.025240 1.1351 2.3333 0.6667 0.5501
## Trifrepe 0.027416 0.017993 1.5237 3.3333 2.1667 0.5968
## Poatriv 0.026724 0.027452 0.9735 4.6667 5.5000 0.6424
## Anthodor 0.026238 0.020621 1.2724 2.1667 0.0000 0.6872
## Juncbufo 0.024559 0.024039 1.0216 0.3333 1.8333 0.7291
## Bromhord 0.024328 0.019807 1.2283 2.0000 0.5000 0.7706
## Bracruta 0.021988 0.021354 1.0297 2.0000 2.0000 0.8081
## Scorautu 0.018939 0.013653 1.3872 2.8333 2.1667 0.8404
## Trifprat 0.017870 0.021259 0.8406 1.5000 0.0000 0.8709
## Poaprat 0.017839 0.019132 0.9324 3.5000 3.1667 0.9013
## Juncarti 0.016978 0.024765 0.6856 0.0000 1.3333 0.9303
## Bellpere 0.014809 0.014581 1.0156 1.1667 0.6667 0.9555
## Ranuflam 0.009039 0.013245 0.6824 0.0000 0.6667 0.9709
## Eleopalu 0.008600 0.019773 0.4349 0.0000 0.6667 0.9856
## Cirsarve 0.004034 0.009258 0.4357 0.0000 0.3333 0.9925
## Chenalbu 0.002369 0.005468 0.4333 0.0000 0.1667 0.9965
## Vicilath 0.002033 0.004618 0.4403 0.1667 0.0000 1.0000
## Airaprae 0.000000 0.000000 NaN 0.0000 0.0000 1.0000
## Comapalu 0.000000 0.000000 NaN 0.0000 0.0000 1.0000
## Empenigr 0.000000 0.000000 NaN 0.0000 0.0000 1.0000
## Hyporadi 0.000000 0.000000 NaN 0.0000 0.0000 1.0000
## Salirepe 0.000000 0.000000 NaN 0.0000 0.0000 1.0000
## Callcusp 0.000000 0.000000 NaN 0.0000 0.0000 1.0000
##
## Contrast: 3_1
##
## average sd ratio ava avb cumsum
## Poatriv 0.071134 0.02222 3.2013 4.6667 0.00 0.1104
## Lolipere 0.065842 0.04953 1.3292 5.3333 2.25 0.2126
## Bracruta 0.042351 0.03280 1.2912 2.0000 3.25 0.2783
## Scorautu 0.037921 0.03176 1.1939 2.8333 4.50 0.3372
## Planlanc 0.036828 0.02270 1.6226 3.0000 2.00 0.3944
## Trifrepe 0.036651 0.02824 1.2979 3.3333 1.75 0.4512
## Hyporadi 0.035220 0.02812 1.2524 0.0000 2.25 0.5059
## Elymrepe 0.035093 0.03850 0.9114 2.0000 0.00 0.5604
## Rumeacet 0.033387 0.03621 0.9221 2.3333 0.00 0.6122
## Anthodor 0.032373 0.03200 1.0118 2.1667 2.00 0.6624
## Poaprat 0.030417 0.02631 1.1560 3.5000 2.00 0.7097
## Bromhord 0.029558 0.02489 1.1873 2.0000 0.00 0.7555
## Achimill 0.028319 0.01536 1.8441 2.3333 0.50 0.7995
## Salirepe 0.022678 0.02398 0.9458 0.0000 1.50 0.8347
## Trifprat 0.021280 0.02555 0.8330 1.5000 0.00 0.8677
## Airaprae 0.020683 0.02216 0.9333 0.0000 1.25 0.8998
## Sagiproc 0.018199 0.01995 0.9123 0.0000 1.25 0.9281
## Bellpere 0.017753 0.01816 0.9775 1.1667 0.50 0.9556
## Vicilath 0.011242 0.01211 0.9287 0.1667 0.75 0.9731
## Empenigr 0.007323 0.01323 0.5536 0.0000 0.50 0.9844
## Juncbufo 0.005090 0.01172 0.4345 0.3333 0.00 0.9923
## Alopgeni 0.004937 0.01136 0.4347 0.3333 0.00 1.0000
## Agrostol 0.000000 0.00000 NaN 0.0000 0.00 1.0000
## Chenalbu 0.000000 0.00000 NaN 0.0000 0.00 1.0000
## Cirsarve 0.000000 0.00000 NaN 0.0000 0.00 1.0000
## Comapalu 0.000000 0.00000 NaN 0.0000 0.00 1.0000
## Eleopalu 0.000000 0.00000 NaN 0.0000 0.00 1.0000
## Juncarti 0.000000 0.00000 NaN 0.0000 0.00 1.0000
## Ranuflam 0.000000 0.00000 NaN 0.0000 0.00 1.0000
## Callcusp 0.000000 0.00000 NaN 0.0000 0.00 1.0000
##
## Contrast: 3_4
##
## average sd ratio ava avb cumsum
## Lolipere 0.084704 0.038558 2.1968 5.3333 0.00 0.09603
## Eleopalu 0.080308 0.026995 2.9749 0.0000 5.25 0.18707
## Agrostol 0.076247 0.020401 3.7374 0.0000 5.00 0.27351
## Poatriv 0.061858 0.026316 2.3506 4.6667 0.50 0.34364
## Poaprat 0.054987 0.019512 2.8181 3.5000 0.00 0.40598
## Trifrepe 0.044678 0.034197 1.3065 3.3333 1.75 0.45663
## Planlanc 0.042303 0.032421 1.3048 3.0000 0.00 0.50458
## Bracruta 0.038932 0.025172 1.5466 2.0000 3.00 0.54872
## Ranuflam 0.038285 0.014364 2.6654 0.0000 2.50 0.59212
## Callcusp 0.038035 0.025820 1.4731 0.0000 2.50 0.63524
## Juncarti 0.037813 0.024529 1.5416 0.0000 2.50 0.67811
## Achimill 0.034414 0.012390 2.7776 2.3333 0.00 0.71712
## Elymrepe 0.033755 0.036292 0.9301 2.0000 0.00 0.75539
## Rumeacet 0.032477 0.034996 0.9280 2.3333 0.00 0.79221
## Anthodor 0.030478 0.024163 1.2614 2.1667 0.00 0.82676
## Bromhord 0.028722 0.023955 1.1990 2.0000 0.00 0.85932
## Scorautu 0.027794 0.017188 1.6170 2.8333 1.50 0.89083
## Trifprat 0.020706 0.024733 0.8372 1.5000 0.00 0.91431
## Salirepe 0.018308 0.033071 0.5536 0.0000 1.25 0.93506
## Alopgeni 0.016792 0.024399 0.6882 0.3333 1.00 0.95410
## Bellpere 0.016654 0.017807 0.9352 1.1667 0.00 0.97298
## Comapalu 0.016530 0.017597 0.9393 0.0000 1.00 0.99172
## Juncbufo 0.004940 0.011311 0.4367 0.3333 0.00 0.99732
## Vicilath 0.002365 0.005413 0.4368 0.1667 0.00 1.00000
## Airaprae 0.000000 0.000000 NaN 0.0000 0.00 1.00000
## Chenalbu 0.000000 0.000000 NaN 0.0000 0.00 1.00000
## Cirsarve 0.000000 0.000000 NaN 0.0000 0.00 1.00000
## Empenigr 0.000000 0.000000 NaN 0.0000 0.00 1.00000
## Hyporadi 0.000000 0.000000 NaN 0.0000 0.00 1.00000
## Sagiproc 0.000000 0.000000 NaN 0.0000 0.00 1.00000
##
## Contrast: 2_1
##
## average sd ratio ava avb cumsum
## Poatriv 0.086332 0.034365 2.5122 5.5000 0.00 0.1160
## Alopgeni 0.079013 0.037490 2.1076 5.0000 0.00 0.2222
## Agrostol 0.071974 0.024236 2.9697 4.6667 0.00 0.3189
## Lolipere 0.045515 0.035550 1.2803 2.8333 2.25 0.3801
## Bracruta 0.036395 0.024788 1.4683 2.0000 3.25 0.4290
## Scorautu 0.035223 0.020888 1.6863 2.1667 4.50 0.4763
## Elymrepe 0.034442 0.037075 0.9290 2.3333 0.00 0.5226
## Hyporadi 0.034208 0.026214 1.3050 0.0000 2.25 0.5686
## Anthodor 0.032868 0.034303 0.9582 0.0000 2.00 0.6127
## Poaprat 0.032448 0.023611 1.3743 3.1667 2.00 0.6563
## Planlanc 0.031236 0.018827 1.6591 0.0000 2.00 0.6983
## Sagiproc 0.030263 0.024672 1.2266 2.5000 1.25 0.7390
## Juncbufo 0.029473 0.030854 0.9552 1.8333 0.00 0.7786
## Salirepe 0.022104 0.022682 0.9745 0.0000 1.50 0.8083
## Juncarti 0.020055 0.029245 0.6858 1.3333 0.00 0.8352
## Airaprae 0.020009 0.020637 0.9696 0.0000 1.25 0.8621
## Trifrepe 0.016639 0.017363 0.9583 2.1667 1.75 0.8845
## Bellpere 0.012598 0.015378 0.8192 0.6667 0.50 0.9014
## Vicilath 0.010842 0.012071 0.8982 0.0000 0.75 0.9160
## Ranuflam 0.010799 0.015811 0.6830 0.6667 0.00 0.9305
## Rumeacet 0.010455 0.015293 0.6836 0.6667 0.00 0.9445
## Eleopalu 0.010180 0.023430 0.4345 0.6667 0.00 0.9582
## Achimill 0.009284 0.016490 0.5630 0.0000 0.50 0.9707
## Empenigr 0.007150 0.012680 0.5639 0.0000 0.50 0.9803
## Bromhord 0.007088 0.016294 0.4350 0.5000 0.00 0.9898
## Cirsarve 0.004725 0.010863 0.4350 0.3333 0.00 0.9962
## Chenalbu 0.002855 0.006584 0.4336 0.1667 0.00 1.0000
## Comapalu 0.000000 0.000000 NaN 0.0000 0.00 1.0000
## Trifprat 0.000000 0.000000 NaN 0.0000 0.00 1.0000
## Callcusp 0.000000 0.000000 NaN 0.0000 0.00 1.0000
##
## Contrast: 2_4
##
## average sd ratio ava avb cumsum
## Poatriv 0.076730 0.035297 2.1738 5.5000 0.50 0.1165
## Eleopalu 0.068316 0.032100 2.1282 0.6667 5.25 0.2202
## Alopgeni 0.065917 0.039379 1.6739 5.0000 1.00 0.3203
## Poaprat 0.046536 0.024779 1.8781 3.1667 0.00 0.3910
## Lolipere 0.040993 0.034867 1.1757 2.8333 0.00 0.4532
## Sagiproc 0.037422 0.023982 1.5605 2.5000 0.00 0.5101
## Callcusp 0.037049 0.023751 1.5599 0.0000 2.50 0.5663
## Trifrepe 0.035339 0.018937 1.8661 2.1667 1.75 0.6200
## Elymrepe 0.033468 0.035748 0.9362 2.3333 0.00 0.6708
## Juncarti 0.032485 0.024290 1.3373 1.3333 2.50 0.7201
## Bracruta 0.030206 0.019581 1.5426 2.0000 3.00 0.7660
## Juncbufo 0.028542 0.029595 0.9644 1.8333 0.00 0.8093
## Ranuflam 0.026823 0.018649 1.4383 0.6667 2.50 0.8501
## Agrostol 0.021487 0.019803 1.0851 4.6667 5.00 0.8827
## Salirepe 0.017875 0.031700 0.5639 0.0000 1.25 0.9098
## Comapalu 0.016027 0.016444 0.9746 0.0000 1.00 0.9342
## Rumeacet 0.010135 0.014719 0.6885 0.6667 0.00 0.9496
## Bellpere 0.009538 0.013832 0.6896 0.6667 0.00 0.9641
## Scorautu 0.009421 0.013442 0.7008 2.1667 1.50 0.9784
## Bromhord 0.006897 0.015787 0.4369 0.5000 0.00 0.9888
## Cirsarve 0.004598 0.010525 0.4369 0.3333 0.00 0.9958
## Chenalbu 0.002757 0.006317 0.4365 0.1667 0.00 1.0000
## Achimill 0.000000 0.000000 NaN 0.0000 0.00 1.0000
## Airaprae 0.000000 0.000000 NaN 0.0000 0.00 1.0000
## Anthodor 0.000000 0.000000 NaN 0.0000 0.00 1.0000
## Empenigr 0.000000 0.000000 NaN 0.0000 0.00 1.0000
## Hyporadi 0.000000 0.000000 NaN 0.0000 0.00 1.0000
## Planlanc 0.000000 0.000000 NaN 0.0000 0.00 1.0000
## Trifprat 0.000000 0.000000 NaN 0.0000 0.00 1.0000
## Vicilath 0.000000 0.000000 NaN 0.0000 0.00 1.0000
##
## Contrast: 1_4
##
## average sd ratio ava avb cumsum
## Eleopalu 0.098265 0.03006 3.2686 0.00 5.25 0.1177
## Agrostol 0.093231 0.02119 4.4001 0.00 5.00 0.2294
## Scorautu 0.052080 0.02919 1.7844 4.50 1.50 0.2917
## Ranuflam 0.046847 0.01632 2.8700 0.00 2.50 0.3479
## Callcusp 0.046477 0.03138 1.4810 0.00 2.50 0.4035
## Juncarti 0.046157 0.02950 1.5647 0.00 2.50 0.4588
## Trifrepe 0.042804 0.03664 1.1683 1.75 1.75 0.5101
## Bracruta 0.042329 0.03965 1.0675 3.25 3.00 0.5608
## Hyporadi 0.041620 0.03195 1.3028 2.25 0.00 0.6106
## Anthodor 0.040746 0.04350 0.9366 2.00 0.00 0.6594
## Lolipere 0.038632 0.04988 0.7744 2.25 0.00 0.7057
## Planlanc 0.038219 0.02329 1.6412 2.00 0.00 0.7515
## Poaprat 0.036516 0.02753 1.3265 2.00 0.00 0.7952
## Salirepe 0.036288 0.03376 1.0749 1.50 1.25 0.8387
## Airaprae 0.024651 0.02570 0.9594 1.25 0.00 0.8682
## Sagiproc 0.021247 0.02299 0.9241 1.25 0.00 0.8936
## Comapalu 0.020538 0.02173 0.9453 0.00 1.00 0.9182
## Alopgeni 0.017127 0.03097 0.5530 0.00 1.00 0.9388
## Vicilath 0.013007 0.01462 0.8898 0.75 0.00 0.9543
## Achimill 0.011816 0.02128 0.5553 0.50 0.00 0.9685
## Bellpere 0.009189 0.01651 0.5567 0.50 0.00 0.9795
## Poatriv 0.008564 0.01549 0.5530 0.00 0.50 0.9898
## Empenigr 0.008557 0.01536 0.5570 0.50 0.00 1.0000
## Bromhord 0.000000 0.00000 NaN 0.00 0.00 1.0000
## Chenalbu 0.000000 0.00000 NaN 0.00 0.00 1.0000
## Cirsarve 0.000000 0.00000 NaN 0.00 0.00 1.0000
## Elymrepe 0.000000 0.00000 NaN 0.00 0.00 1.0000
## Juncbufo 0.000000 0.00000 NaN 0.00 0.00 1.0000
## Rumeacet 0.000000 0.00000 NaN 0.00 0.00 1.0000
## Trifprat 0.000000 0.00000 NaN 0.00 0.00 1.0000
## Permutation: free
## Number of permutations: 0sim2 <- simper(dune, dune.env.original$Management) #try management groupings
summary(sim2)##
## Contrast: SF_BF
##
## average sd ratio ava avb cumsum
## Agrostol 0.061374 0.034193 1.7949 4.6667 0.0000 0.09824
## Alopgeni 0.052667 0.036476 1.4439 4.3333 0.6667 0.18255
## Lolipere 0.048116 0.039445 1.2198 3.0000 6.0000 0.25957
## Trifrepe 0.046297 0.025525 1.8138 1.3333 4.6667 0.33368
## Poatriv 0.046020 0.033801 1.3615 4.6667 3.6667 0.40734
## Scorautu 0.043697 0.024922 1.7534 1.3333 4.3333 0.47729
## Bromhord 0.033677 0.025860 1.3023 0.5000 2.6667 0.53120
## Achimill 0.030152 0.020821 1.4482 0.1667 2.3333 0.57947
## Planlanc 0.028585 0.021549 1.3265 0.0000 2.0000 0.62522
## Elymrepe 0.028074 0.029778 0.9428 2.0000 1.3333 0.67016
## Bracruta 0.025501 0.023902 1.0669 2.0000 2.0000 0.71098
## Poaprat 0.025129 0.023967 1.0485 2.5000 4.0000 0.75121
## Sagiproc 0.024326 0.022149 1.0983 1.8333 0.6667 0.79014
## Bellpere 0.019859 0.017088 1.1622 0.6667 1.6667 0.82193
## Eleopalu 0.018611 0.042958 0.4333 1.3333 0.0000 0.85172
## Anthodor 0.017543 0.025804 0.6798 0.0000 1.3333 0.87981
## Juncbufo 0.016031 0.023708 0.6762 1.1667 0.0000 0.90547
## Vicilath 0.014671 0.013306 1.1026 0.0000 1.0000 0.92895
## Hyporadi 0.010286 0.015198 0.6768 0.0000 0.6667 0.94542
## Ranuflam 0.009306 0.013595 0.6845 0.6667 0.0000 0.96031
## Juncarti 0.006979 0.016109 0.4333 0.5000 0.0000 0.97148
## Callcusp 0.006979 0.016109 0.4333 0.5000 0.0000 0.98266
## Rumeacet 0.004526 0.010444 0.4333 0.3333 0.0000 0.98990
## Cirsarve 0.003983 0.009185 0.4336 0.3333 0.0000 0.99628
## Chenalbu 0.002326 0.005370 0.4333 0.1667 0.0000 1.00000
## Airaprae 0.000000 0.000000 NaN 0.0000 0.0000 1.00000
## Comapalu 0.000000 0.000000 NaN 0.0000 0.0000 1.00000
## Empenigr 0.000000 0.000000 NaN 0.0000 0.0000 1.00000
## Salirepe 0.000000 0.000000 NaN 0.0000 0.0000 1.00000
## Trifprat 0.000000 0.000000 NaN 0.0000 0.0000 1.00000
##
## Contrast: SF_HF
##
## average sd ratio ava avb cumsum
## Agrostol 0.047380 0.031273 1.5151 4.6667 1.4 0.08351
## Alopgeni 0.046433 0.032897 1.4115 4.3333 1.6 0.16535
## Lolipere 0.041986 0.027007 1.5546 3.0000 4.0 0.23935
## Planlanc 0.039198 0.033208 1.1804 0.0000 3.0 0.30844
## Rumeacet 0.038992 0.027369 1.4247 0.3333 3.2 0.37716
## Elymrepe 0.031877 0.029550 1.0787 2.0000 2.0 0.43334
## Poatriv 0.028466 0.021522 1.3227 4.6667 4.8 0.48352
## Bracruta 0.025261 0.021044 1.2004 2.0000 2.8 0.52804
## Eleopalu 0.024974 0.038877 0.6424 1.3333 0.8 0.57206
## Poaprat 0.023932 0.019180 1.2478 2.5000 3.4 0.61424
## Anthodor 0.023409 0.021430 1.0923 0.0000 1.8 0.65550
## Sagiproc 0.023144 0.020479 1.1301 1.8333 0.8 0.69629
## Trifprat 0.023080 0.023432 0.9850 0.0000 1.8 0.73697
## Juncarti 0.022850 0.025677 0.8899 0.5000 1.6 0.77724
## Trifrepe 0.022383 0.019487 1.1486 1.3333 2.8 0.81669
## Juncbufo 0.021643 0.022237 0.9733 1.1667 1.2 0.85484
## Scorautu 0.020509 0.016422 1.2489 1.3333 2.8 0.89099
## Achimill 0.015183 0.011393 1.3326 0.1667 1.2 0.91775
## Bromhord 0.013375 0.014504 0.9222 0.5000 0.8 0.94132
## Ranuflam 0.010661 0.013387 0.7964 0.6667 0.4 0.96011
## Bellpere 0.009991 0.012571 0.7948 0.6667 0.4 0.97772
## Callcusp 0.006623 0.015076 0.4393 0.5000 0.0 0.98939
## Cirsarve 0.003809 0.008669 0.4394 0.3333 0.0 0.99611
## Chenalbu 0.002208 0.005025 0.4393 0.1667 0.0 1.00000
## Airaprae 0.000000 0.000000 NaN 0.0000 0.0 1.00000
## Comapalu 0.000000 0.000000 NaN 0.0000 0.0 1.00000
## Empenigr 0.000000 0.000000 NaN 0.0000 0.0 1.00000
## Hyporadi 0.000000 0.000000 NaN 0.0000 0.0 1.00000
## Salirepe 0.000000 0.000000 NaN 0.0000 0.0 1.00000
## Vicilath 0.000000 0.000000 NaN 0.0000 0.0 1.00000
##
## Contrast: SF_NM
##
## average sd ratio ava avb cumsum
## Poatriv 0.078284 0.040947 1.9118 4.6667 0.0000 0.1014
## Alopgeni 0.071219 0.046958 1.5167 4.3333 0.0000 0.1936
## Agrostol 0.056508 0.044176 1.2792 4.6667 2.1667 0.2667
## Lolipere 0.054851 0.059914 0.9155 3.0000 0.3333 0.3378
## Eleopalu 0.048027 0.047168 1.0182 1.3333 2.1667 0.3999
## Poaprat 0.040724 0.031790 1.2810 2.5000 0.6667 0.4527
## Bracruta 0.040008 0.034398 1.1631 2.0000 2.8333 0.5045
## Elymrepe 0.035598 0.038515 0.9243 2.0000 0.0000 0.5506
## Scorautu 0.032475 0.034813 0.9328 1.3333 3.1667 0.5926
## Trifrepe 0.030430 0.031634 0.9619 1.3333 1.8333 0.6320
## Sagiproc 0.030304 0.030477 0.9943 1.8333 0.5000 0.6712
## Salirepe 0.029275 0.032014 0.9144 0.0000 1.8333 0.7092
## Anthodor 0.024541 0.036694 0.6688 0.0000 1.3333 0.7409
## Callcusp 0.022763 0.029443 0.7731 0.5000 1.1667 0.7704
## Ranuflam 0.022566 0.022819 0.9889 0.6667 1.3333 0.7996
## Juncarti 0.022543 0.028598 0.7883 0.5000 1.1667 0.8288
## Hyporadi 0.020108 0.031291 0.6426 0.0000 1.1667 0.8548
## Juncbufo 0.019860 0.029034 0.6840 1.1667 0.0000 0.8806
## Planlanc 0.015420 0.022769 0.6772 0.0000 0.8333 0.9005
## Airaprae 0.014883 0.021881 0.6802 0.0000 0.8333 0.9198
## Bellpere 0.012317 0.015921 0.7737 0.6667 0.3333 0.9357
## Comapalu 0.011883 0.017407 0.6826 0.0000 0.6667 0.9511
## Achimill 0.009294 0.014931 0.6224 0.1667 0.3333 0.9632
## Bromhord 0.007172 0.016333 0.4391 0.5000 0.0000 0.9724
## Rumeacet 0.005590 0.012751 0.4384 0.3333 0.0000 0.9797
## Empenigr 0.005225 0.012001 0.4354 0.0000 0.3333 0.9864
## Cirsarve 0.004782 0.010889 0.4391 0.3333 0.0000 0.9926
## Chenalbu 0.002893 0.006602 0.4382 0.1667 0.0000 0.9964
## Vicilath 0.002792 0.006425 0.4345 0.0000 0.1667 1.0000
## Trifprat 0.000000 0.000000 NaN 0.0000 0.0000 1.0000
##
## Contrast: BF_HF
##
## average sd ratio ava avb cumsum
## Rumeacet 0.038666 0.02606 1.4838 0.0000 3.2 0.08163
## Poatriv 0.033301 0.02579 1.2911 3.6667 4.8 0.15194
## Planlanc 0.031852 0.01830 1.7401 2.0000 3.0 0.21918
## Bromhord 0.028651 0.01799 1.5926 2.6667 0.8 0.27967
## Lolipere 0.028431 0.02215 1.2834 6.0000 4.0 0.33970
## Elymrepe 0.027822 0.02959 0.9404 1.3333 2.0 0.39843
## Trifrepe 0.025838 0.01656 1.5603 4.6667 2.8 0.45298
## Anthodor 0.023582 0.02042 1.1547 1.3333 1.8 0.50277
## Achimill 0.023426 0.01474 1.5893 2.3333 1.2 0.55223
## Bracruta 0.022733 0.01802 1.2617 2.0000 2.8 0.60022
## Alopgeni 0.021610 0.02308 0.9363 0.6667 1.6 0.64584
## Trifprat 0.021514 0.02207 0.9747 0.0000 1.8 0.69126
## Juncarti 0.020084 0.02555 0.7860 0.0000 1.6 0.73367
## Scorautu 0.019318 0.01356 1.4241 4.3333 2.8 0.77445
## Bellpere 0.018290 0.01486 1.2305 1.6667 0.4 0.81306
## Agrostol 0.017605 0.02284 0.7708 0.0000 1.4 0.85023
## Juncbufo 0.015000 0.02066 0.7260 0.0000 1.2 0.88190
## Vicilath 0.012848 0.01140 1.1274 1.0000 0.0 0.90903
## Sagiproc 0.011685 0.01297 0.9008 0.6667 0.8 0.93369
## Eleopalu 0.010169 0.02111 0.4817 0.0000 0.8 0.95516
## Hyporadi 0.008950 0.01312 0.6824 0.6667 0.0 0.97406
## Poaprat 0.007203 0.01010 0.7133 4.0000 3.4 0.98927
## Ranuflam 0.005084 0.01055 0.4817 0.0000 0.4 1.00000
## Airaprae 0.000000 0.00000 NaN 0.0000 0.0 1.00000
## Chenalbu 0.000000 0.00000 NaN 0.0000 0.0 1.00000
## Cirsarve 0.000000 0.00000 NaN 0.0000 0.0 1.00000
## Comapalu 0.000000 0.00000 NaN 0.0000 0.0 1.00000
## Empenigr 0.000000 0.00000 NaN 0.0000 0.0 1.00000
## Salirepe 0.000000 0.00000 NaN 0.0000 0.0 1.00000
## Callcusp 0.000000 0.00000 NaN 0.0000 0.0 1.00000
##
## Contrast: BF_NM
##
## average sd ratio ava avb cumsum
## Lolipere 0.090681 0.02644 3.4292 6.0000 0.3333 0.1243
## Poatriv 0.054684 0.04465 1.2248 3.6667 0.0000 0.1992
## Poaprat 0.052511 0.01813 2.8966 4.0000 0.6667 0.2712
## Trifrepe 0.051287 0.02756 1.8611 4.6667 1.8333 0.3415
## Bromhord 0.039689 0.02920 1.3590 2.6667 0.0000 0.3959
## Bracruta 0.035723 0.02869 1.2452 2.0000 2.8333 0.4448
## Eleopalu 0.033759 0.03573 0.9449 0.0000 2.1667 0.4911
## Agrostol 0.033446 0.03473 0.9630 0.0000 2.1667 0.5369
## Achimill 0.033190 0.02338 1.4198 2.3333 0.3333 0.5824
## Scorautu 0.031356 0.02026 1.5480 4.3333 3.1667 0.6254
## Anthodor 0.028060 0.03295 0.8517 1.3333 1.3333 0.6638
## Planlanc 0.027319 0.02193 1.2458 2.0000 0.8333 0.7013
## Salirepe 0.026770 0.02927 0.9145 0.0000 1.8333 0.7379
## Bellpere 0.023529 0.01909 1.2322 1.6667 0.3333 0.7702
## Hyporadi 0.021721 0.02450 0.8864 0.6667 1.1667 0.8000
## Ranuflam 0.020314 0.02275 0.8928 0.0000 1.3333 0.8278
## Elymrepe 0.019993 0.02926 0.6833 1.3333 0.0000 0.8552
## Callcusp 0.017833 0.02681 0.6653 0.0000 1.1667 0.8796
## Juncarti 0.017694 0.02600 0.6806 0.0000 1.1667 0.9039
## Vicilath 0.015773 0.01447 1.0902 1.0000 0.1667 0.9255
## Sagiproc 0.015432 0.01857 0.8310 0.6667 0.5000 0.9466
## Airaprae 0.013410 0.01969 0.6809 0.0000 0.8333 0.9650
## Comapalu 0.010739 0.01571 0.6835 0.0000 0.6667 0.9797
## Alopgeni 0.009997 0.01463 0.6833 0.6667 0.0000 0.9934
## Empenigr 0.004787 0.01105 0.4332 0.0000 0.3333 1.0000
## Chenalbu 0.000000 0.00000 NaN 0.0000 0.0000 1.0000
## Cirsarve 0.000000 0.00000 NaN 0.0000 0.0000 1.0000
## Juncbufo 0.000000 0.00000 NaN 0.0000 0.0000 1.0000
## Rumeacet 0.000000 0.00000 NaN 0.0000 0.0000 1.0000
## Trifprat 0.000000 0.00000 NaN 0.0000 0.0000 1.0000
##
## Contrast: HF_NM
##
## average sd ratio ava avb cumsum
## Poatriv 0.071553 0.013681 5.2302 4.8 0.0000 0.09913
## Lolipere 0.054533 0.029625 1.8408 4.0 0.3333 0.17468
## Rumeacet 0.046546 0.030806 1.5109 3.2 0.0000 0.23917
## Poaprat 0.041750 0.018852 2.2146 3.4 0.6667 0.29701
## Planlanc 0.041633 0.029560 1.4084 3.0 0.8333 0.35469
## Bracruta 0.035340 0.020104 1.7579 2.8 2.8333 0.40365
## Eleopalu 0.032043 0.032315 0.9916 0.8 2.1667 0.44805
## Agrostol 0.031915 0.028887 1.1048 1.4 2.1667 0.49227
## Trifrepe 0.030372 0.022871 1.3280 2.8 1.8333 0.53434
## Elymrepe 0.029811 0.038676 0.7708 2.0 0.0000 0.57565
## Anthodor 0.028717 0.024799 1.1580 1.8 1.3333 0.61543
## Juncarti 0.027414 0.028537 0.9607 1.6 1.1667 0.65341
## Trifprat 0.025843 0.025972 0.9951 1.8 0.0000 0.68922
## Salirepe 0.025340 0.027291 0.9285 0.0 1.8333 0.72432
## Alopgeni 0.024459 0.032399 0.7549 1.6 0.0000 0.75821
## Scorautu 0.020703 0.014125 1.4658 2.8 3.1667 0.78689
## Ranuflam 0.019285 0.019939 0.9672 0.4 1.3333 0.81361
## Juncbufo 0.018181 0.024648 0.7376 1.2 0.0000 0.83880
## Hyporadi 0.017141 0.026548 0.6457 0.0 1.1667 0.86255
## Callcusp 0.016833 0.024901 0.6760 0.0 1.1667 0.88587
## Achimill 0.016555 0.014900 1.1111 1.2 0.3333 0.90881
## Sagiproc 0.015282 0.016535 0.9243 0.8 0.5000 0.92998
## Airaprae 0.012605 0.018243 0.6910 0.0 0.8333 0.94744
## Bromhord 0.012094 0.015169 0.7973 0.8 0.0000 0.96420
## Comapalu 0.010105 0.014556 0.6942 0.0 0.6667 0.97820
## Bellpere 0.008801 0.013732 0.6409 0.4 0.3333 0.99039
## Empenigr 0.004536 0.010325 0.4393 0.0 0.3333 0.99668
## Vicilath 0.002399 0.005461 0.4393 0.0 0.1667 1.00000
## Chenalbu 0.000000 0.000000 NaN 0.0 0.0000 1.00000
## Cirsarve 0.000000 0.000000 NaN 0.0 0.0000 1.00000
## Permutation: free
## Number of permutations: 0
##Q29
Question 29: Diagnostic species for a priori defined groups
Use the resulting data table to find the two top indicator species for each of the four Management types!
In this exercise we will investigate if any particular species are indicative for the four different management types. The CANOCO version of this analysis is a simplified version of these types of tests. The full versions also give significance testing and more comprehensive result presentations. See for instance IndVal by Dufresne & Legendre. 1. Go to menu Data /Add new table(s) / Indicator values. 2. Mark the Species data table, and Management as the factor defining groups. 3. Opt for Indicator Value (quantitative). 4. Export the table to Excel to be able to sort the data to find the highest values = best indicators for a group.
library(labdsv)
library(tibble)
iva <- indval(dune, dune.env.original$Management) #this time we don't need dummy variables for management
iva.df <- as.data.frame(iva$indval)
#arrange in descending order to find best indicator species, change BF to other management types as needed
arrange(rownames_to_column(iva.df), desc(HF))##Q30
Question 30: Molecular Data
Molecular
#to follow, Stefan
##Q31
Question 31: PLS
Examine the p-values to explore which response variables actually give a significant result.
Unless you are already familiar with PLS in R we suggest that you follow the exercise in SIMCA and then try performing the same analysis by modifying the R code below. Remember you can always read the helpfile for any R function by entering ? and the function name in the console, such as ?plsr in this case.
data.pls <- readRDS("Trend_Lakes_2015_PLS.RDS")
#View(data.pls)
library(pls)
#PLS regression, like PCA, seeks to find components which maximize the variability of predictors but differs from PCA as PLS requires #the components to have maximum correlation with the response.
#Change 'Antal arter' (number of species in Swedish) to your choice of response variable(s)
pls.fit <- plsr(`Antal arter` ~ ., data = data.pls[c(9:35)], scale = TRUE, validation = "CV")
#The predictor variables are mapped to a smaller set of variables, and within that smaller space we perform a regression against the #outcome variable. PLS aims to choose new mapped variables that maximally explain the outcome variable.
#The next step is to remove unwanted variables and then build a model.
#Cross validation is used to find the optimal number of retained dimensions.
#Then the model is rebuilt with this optimal number of dimensions.
#Find the number of dimensions with lowest cross validation error
cv <- RMSEP(pls.fit)
best.dims <- which.min(cv$val[estimate = "adjCV", , ]) - 1
best.dims ## 2 comps
## 2# Rerun the model
pls.fit2 <- plsr(`Antal arter` ~ ., data = data.pls[c(6:35)], ncomp = best.dims)
#Finally, we extract the useful information and format the output.
coefficients <- coef(pls.fit2)
sum.coef <- sum(sapply(coefficients, abs))
coefficients <- coefficients * 100 / sum.coef
coefficients <- sort(coefficients[, 1 , 1])
barplot(tail(coefficients, 5))
#The regression coefficients are normalized so their absolute sum is 100 and the result is sorted.
#The results below show that 'Tot-N_TNb' and 'Mn' are the two most importnat positive predictors of 'Antal arter'. You could run the code #barplot(head(coefficients, 5)) to see that at the other end of the scale for negative predictors.
summary(pls.fit2)## Data: X dimension: 32 29
## Y dimension: 32 1
## Fit method: kernelpls
## Number of components considered: 2
## TRAINING: % variance explained
## 1 comps 2 comps
## X 61.778 98.653
## Antal arter 7.295 8.371plot(pls.fit)
validationplot(pls.fit2, val.type = "MSEP")
predplot(pls.fit2)
coefplot(pls.fit2)